Featured Post

#64 The oscillation in peripheral blood-II 血球的數量消長_中



In the previous episode, we have discussed the basic knowledge about “cyclic neutropenia.” Today we will use MATLAB to simulate the dynamical behavior observed in cyclic neutropenia.

在上集我們已經簡單介紹過「週期性嗜中性球低下」的基本知識(請見臉書粉絲專頁),這集我們要來用程式模擬重現出「週期性嗜中性球低下」的現象。



Let’s briefly review the hematopoiesis of neutrophils (We actually skipped some steps.):
先複習一下嗜中性球的製造過程:(事實上還是跳過了一點點步驟)
Hematopoietic stem cell (HSC) → Myeloblast (MyB)→ Promyelocyte (ProM)→ Myelocyte (MyC)→ Metamyelocyte (MetaM) → Neutrophil (NO)。
翻譯成中文就是:
造血幹細胞→骨髓芽細胞→前骨髓細胞→骨髓細胞→後骨髓細胞→嗜中性球。
造血幹細胞要往骨髓芽細胞分化需要G-CSF的刺激。




Clinically the physicians discovered that the peripheral blood of patients suffered from cyclic neutropenia possess the following properties:
(1) The number of neutrophils varies with a fixed time period.
(2) G-CSF also varies with a fixed time period. However, the oscillation of neutrophil and G-CSF are in antiphase.
(3) The number of erythrocytes and platelets also vary with a fixed time period. The time period is comparable to the time period of neutrophil oscillation.
(4) The oscillation of erythrocytes & platelets and the oscillation of neutrophil are in antiphase.
(5) The erythropoietin varies with a fixed time period.
(6) The oscillation of erythropoietin and erythrocytes are in-phase.
在臨床上前輩觀察到「週期性嗜中性球低下」患者的周邊血液有以下現象:
(1) 嗜中性球的數目有週期性變化,週期固定。
(2)  G-CSF也有週期性變化,但G-CSF的週期性變化與嗜中性球的週期性變化是anti-phase。
(3) 患者的紅血球製造與血小板製造也有週期性變化,且週期和嗜中性球的週期相若。
(4)  患者的紅血球製造與血小板製造的週期性變化與嗜中性球的週期性變化也是anti-phase。
(5) 患者的EPO(紅血球生成素,刺激造血幹細胞往紅血球分化)也有週期性變化。
(6) 患者的EPO和紅血球的週期性變化是in-phase。



How could these happen? To explain it, John Milton & Michael C Mackey formulated a model:
那究竟是如何產生這些變化的呢?為此前John Milton & Michael C Mackey提出了下面的模型:
 Figure originally from:Milton, J. G. & Mackey, M. C. (1989). Periodic haematological diseases: mystical entities or dynamical disorders? J R Coll Physicians Lond 23(4): 236-241. Figure 1.


This is comparable to our current knowledge about hematopoiesis. The only difference is that the production of G-CSF is regulated by the number of neutrophils in peripheral blood. Is that sufficient to reproduce the oscillation observed in cyclic neutropenia? Let’s do it by ourselves !
其實和我們現今的認知是差不多的,只是G-CSF的製造會受到周邊嗜中性球數量的調控。這樣就足以解釋「週期性嗜中性球低下」的現象了嗎?我們就自己實際來操作一下就知道囉!

Simulation by MATLAB程式模擬


We define a function which stores our system of ODE:
這裡用MATLAB進行模擬,先寫好一個定義微分方程組的函數:



function dx = CyclicNeutropenia(t, x)

      %Simulation of cyclic neutropenia
      %HSC = hematopoietic stem cell
      %MyB = Myeloblast
      %ProM = Promyelocyte, it possess nonspecific granules
      %MyC = Myelocyte, it possess specific granules
      %MetaM = Metamyelocyte, it is post-mitotic
      %NO = Neutrophil
      %GCSF= G-CSF
      dx = zeros(7,1);
      a = 1; f = 0.3; b1 = 0.01; d_HSC = 0.5;
      b2 = 0.7; c1 = 0.5; d_MyB = 1;
      b3 = 0.5; c2 = 0.3; d_ProM = 1;
      b4 = 0.5; c3 = 0.3; d_MyC = 1;
      c4 = 0.3; d_MetaM = 1;
      d_NO = 10;
      GCSF_prod = 100; K = 7500; d_GCSF = 10;
     
      dx(1) = a*(1-f)*x(7)*x(1) + b1*x(1) - d_HSC*x(1);  
      %HSC' = +a*(1-f)*GCSF*HSC + b1*HSC - d_HSC*HSC
      dx(2) = a*f*x(7)*x(1) + b2*x(2) -c1*x(2) - d_MyB*x(2);  
      %MyB' = +a*f*GCSF*HSC + b2*MyB - c1*MyB - d_MyB*MyB
      dx(3) = c1*x(2) + b3*x(3) - c2*x(3) - d_ProM*x(3);   
      %ProM' = +c1*MyB +b3*ProM -c2*ProM - d_ProM*ProM
      dx(4) = c2*x(3) + b4*x(4) - c3*x(4) - d_MyC*x(4); 
      %MyC' = +c2*ProM + b4*MyC - c3*MyC - d_MyC*MyC
      dx(5) = c3*x(4) - c4*x(5) - d_MetaM*x(5);
      %MetaM' = +c3*MyC - c4*MetaM - d_MetaM*MetaM
      dx(6) = c4*x(5) - d_NO*x(6);  
      %NO' = +c4*MetaM - d_NO*NO
      dx(7) = GCSF_prod*(1 - x(6)/(K + x(6))) - d_GCSF*x(7); 
      %GCSF' = +GCSF_prod*(1 - NO/(K + NO)) - d_GCSF*GCSF

end
 



We should always check whether the system of ODE is reasonable. “a” means the responsiveness of HSC to G-CSF. “f” means the fraction of HSC that will differentiate upon the stimulation of G-CSF. “b” means the self-renewability of cells. “d” means the apoptotic rate of cells. The production of G-CSF is regulated by neutrophils. It seems reasonable to me. So let’s run it! We solve the above system with function “ode45.”

 在看參數列表或實際RUN之前,我們要先確認一下我們定義的方程組是不是合理的。a 表示造血幹細胞對G-CSF的反應程度,f 則表示有多少比例會因此刺激而開始往骨髓細胞分化。b 表示細胞有自我複製的能力。d則是這些細胞的死亡速率。G-CSF的製造受到嗜中性球數目的調控。看起來應該是蠻合理的,所以我們來跑跑看吧,用ode45來解的結果是:



Aside from the initial damped oscillation, the variables reached a fixed point and remained there for a long time. This simulates the peripheral blood conditions of most people, in which the relative blood components remain stable for a long time.
除了一開始的小小震盪之外,達到平衡之後就維持在穩定的數值。這和大部分人一樣,雖然血球數目還是會有多多少少的變動,但基本上是在一個穩定的數值附近(注意這裡y座標是log scale喔,所以3倍以內的變化基本上是看不太出來的。)



How could we reproduce the oscillatory behavior Let’s look back to our 7-dimensional ODE system:
那到底要怎樣才會有週期性變化的發生呢??我們回頭看看我們的七個聯立微分方程組:

%HSC' = +a*(1-f)*GCSF*HSC + b1*HSC - d_HSC*HSC
%MyB' = +a*f*GCSF*HSC + b2*MyB - c1*MyB - d_MyB*MyB
%ProM' = +c1*MyB +b3*ProM -c2*ProM - d_ProM*ProM
%MyC' = +c2*ProM + b4*MyC - c3*MyC - d_MyC*MyC  
%MetaM' = +c3*MyC - c4*MetaM - d_MetaM*MetaM 
%NO' = +c4*MetaM - d_NO*NO 
%GCSF' = +GCSF_prod*(1 - NO/(K + NO)) - d_GCSF*GCSF

 

It’s obvious that the ODE of MyB, ProM, MyC, MetaM are simply transducing the dynamical behavior of HSC to NO. In steady state, their values are directly proportional to the number of cells in previous differentiation states. (For example, [ProM]ss = c1/(d_MyC + c3 - b4)*[MyB]ss.) So simply adjusting parameters such as c1, b3, d_ProM would not change the overall dynamical behavior of the system. The key would lie in the ODE involve HSC, NO, and G-CSF. So possible candidates for adjustment are a, f, b1, d_HSC, d_NO, GCSF_prod, K, d_GCSF.
我們可以看得出來,MyB, ProM, MyC, MetaM的微分方程組,基本上只是在把HSC的動態變化傳遞到NO而已。在穩定狀態下(對時間微分為零),他們的穩定數值都和分化的前一步的細胞的穩定數值有直接的正比關係(例如:[ProM]ss = c1/(d_MyC + c3 - b4)*[MyB]ss。)所以調整這些有關的參數例如c1, b3, d_ProM之類的,只會讓最後的穩定數值增加或減少,並不會影響他的動態行為。所以關鍵應該是在HSC, NO, G-CSF的動態變化才對。所以我們可能會去調整的參數只有a, f, b1, d_HSC, d_NO, GCSF_prod, K, d_GCSF這些而已。




The readers are invited to play with the code above or rewrite it yourself in Python, for example. How could we reproduce the oscillation? In fact only the adjustment of a, f, b1, d_HSC would be productive, and we will explain why in the next episode. To say it in biological words, the genetic mutation causing cyclic neutropenia should affect at least one of the followings:
(1)  The responsiveness of HSC to G-CSF.
(2)  The fraction of HSC which differentiate into myeloblast upon the stimulation of G-CSF.
(3)  The renewal rate of HSC.
The apoptotic rate of HSC.
各位可以用這裡提供的程式碼玩玩看,究竟怎樣調整才能夠製造出週期性的變動呢?事實上我們是可以用數學方法知道,只有a, f, b1, d_HSC四個參數,能夠用些微的調整就讓系統從沒有週期震盪變成有週期震盪,但這我們留到下集再討論。翻譯成生物學的語言就是,要造成「週期性嗜中性球低下」的基因突變,必須影響到以下四個中的其中至少一個:
(1) 造血幹細胞對G-CSF的影響程度。
(2) 造血幹細胞受G-CSF刺激後分化成骨髓前細胞的比例。
(3) 造血幹細胞自我更新的速率。
(4) 造血幹細胞死亡的速率。



In fact, the most common mutation causing cyclic neutropenia involves abnormally increased apoptotic rate of hematopoietic stem cells. So in the following we are only going to demonstrate the effect of adjusting d_HSC.
實際生物實驗發現的是, 「週期性嗜中性球低下」患者的造血幹細胞死亡的速率異常增加,所以我們這裡只演示增加造血幹細胞的死亡速率造成的結果:



We adjust the parameters as the following:
首先我們調整參數如下:
      a = 1; f = 0.3; b1 = 0.01; d_HSC = 1;
      b2 = 0.7; c1 = 0.5; d_MyB = 1;
      b3 = 0.5; c2 = 0.3; d_ProM = 1;
      b4 = 0.5; c3 = 0.3; d_MyC = 1;
      c4 = 0.3; d_MetaM = 1;
      d_NO = 10;
      GCSF_prod = 100; K = 7500; d_GCSF = 10;

And the results look like this:
結果:

The number of neutrophil is now oscillating and is in antiphase with the oscillation of G-CSF.
我們發現周邊血液的嗜中性球數量真的出現週期變化,而且和G-CSF的週期變化是anti-phase。但這名患者的嗜中性球數目並不會降到很低的數值。



If we keep increasing the apoptotic rate of HSC:  
如果我們繼續調大造血幹細胞的死亡速率:
      a = 1; f = 0.3; b1 = 0.01; d_HSC = 5;
      b2 = 0.7; c1 = 0.5; d_MyB = 1;
      b3 = 0.5; c2 = 0.3; d_ProM = 1;
      b4 = 0.5; c3 = 0.3; d_MyC = 1;
      c4 = 0.3; d_MetaM = 1;
      d_NO = 10;
      GCSF_prod = 100; K = 7500; d_GCSF = 10;


We’ll get:
我們就會得到:

And now we have a patient with cyclic neutropenia. The state of having his or her peripheral neutrophil count less than <500/μL is called agranulocytosis and it is particularly dangerous.
注意圖是log scale,所以現在嗜中性球的數目已經會震盪到非常接近零的數目了。當嗜中性球數目<500/μL時我們稱為agranulocytosis,處於這種狀態的病人是非常危險的,必須待在正壓隔離病房裡面才行。基本上我們已經重現了 「週期性嗜中性球低下」的基本動態。關於中性球製造與紅血球製造如何couple在一起,我們留給讀者自己修改我們的模型囉。





It’s worth noting that further adjustment could eliminate the oscillation again:

值得注意的是,如果繼續調整,週期性震盪的現象又會消失:

      a = 0.91; f = 0.3; b1 = 0.01; d_HSC = 6;
      b2 = 0.7; c1 = 0.5; d_MyB = 1;
      b3 = 0.5; c2 = 0.3; d_ProM = 1;
      b4 = 0.5; c3 = 0.3; d_MyC = 1;
      c4 = 0.3; d_MetaM = 1;
      d_NO = 10;
      GCSF_prod = 100; K = 7500; d_GCSF = 10;


These patients become always neutropenic.
這種患者的嗜中性球就一直處於很低的狀態。



And let’s call it a day. In the next episode we are going to simplify our model into a 2 dimensional system of ODE to facilitate discussion upon the dynamical behavior of the system.
今天的講解就到這邊,下集我們會再用更簡化的二維模型探討為什麼只有剛剛列的四個參數會影響系統的動態行為。

/*------------Divider------------*/
Reference
Milton, J. G. & Mackey, M. C. (1989). Periodic haematological diseases: mystical entities or dynamical disorders? J R Coll Physicians Lond 23(4): 236-241.

Comments