Systems biology and engineering, SBE

임영일 CAPEC, Department of Chemical Engineering, Technical University of Denmark, DK-2800 Lyngby, Denmark Phone: +45 4525 2804, Fax: +45 4593 2906, Email: lim@kt.dtu.dk 5. 세포군집모델링 (cell population modeling) Systems biology and engineering (SBE) 3 장에서설명한세포모사는단일세포에대한세포내대사 / 조절 / 신호전달경로와반응에관하여다룬다. 하지만, 생물반응기내에는비균일상태의수백만개세포들로이루어져있다. 따라서세포내물질들 (metabolites, DNA, proteins etc.) 의농도등은단일값이라기보다는분포식으로표현되어야한다. 이러한비균일상태의미생물에대한동적거동의모델링을통한생물반응기 (fermentation reactor) 의설계, 모사, 최적화, 제어등을위하여, 세포입자수지식 (cell population balance equation, PBE) 이나세포앙상블모델 (cell ensemble model) 등이유도된다 (Henson, 2003). PBE 모델은주로효모 (S. cervisiae) 의유산소성장과글루코즈제한환경에서 cell-cycle dependent 진동현상을설명하고, cell ensemble model 은 glycolytic 진동현상을설명하기위하여주로이용되어왔다 (Mhaskar and Hjortso, 2002). 본장에서는단일세포의동적거동 (single cell models and behaviors) 과이들세포들의상호작용및분포에관한 cell population dynamics 을이용한생물반응기의모델링과모사에관하여설명한다. Keywords: cell population dynamics, population balance equation, cell ensemble model 5.1 세포군집수지식 (Cell population balance equations) 미생물배양에있어서세포내또는세포외부의특정물질농도 (e.g., cell mass, glucose/ethanol/oxygen concentrations, culture ph, CO 2 ) 에대한자체적인진동이관찰되어진다. 이러한진동현상은생물반응기의안정적인운전과생산성에큰영향을미치며, 이러한진동현상을설명할수있는동적모델이필요하다. 이들동적모델은크게세가지로구분되어진다 (Zhang et al., 2002): i) a structured/unsegregated modeling on the basis of cell metabolite simulation (Jones and Kompala, 1999), ii) a segregated/unstructured modeling based on the cell population balance equation (PBE) coupled to the mass balance of the rate limiting substrate (Zhu et al., 2000; Zhang et al., 2002), and iii) segregated/structured modeling based on PBE coupled to metabolic reactions relevant to extracellular environment (Mhaskar et al., 2002). 첫번째모델 (structured/unsegregated) 은세포대사반응의직접적인결과가진동현상의원인으로설명하고, 두번째모델 (segregated/unstructured) 에서는진동현상은세포대사반응에기인하기보다는수백만세포군집들의동적거동에기인한다. 세번째모델 (segregated/structured) 은위 2 모델의장점을취합한것으로, 세포내대사반응과세포입자분포식을포함하고있다 (Zhang et al., 2002). Structured/unsegregated modeling ( 체계적단일세포모델링 ) 단일세포내 (unsegregated) 복잡한대사 / 조절 / 신호전달경로및반응을체계적으로구축하여 (structured) 세포내부 / 외부물질등의시간에따른진동현상을설명하려는접근방식이다. 이는본보고서의 3장에서설명한세포의동적모델링에상응하는것으로, 매우현학적이나여러상태에서존재하는수백만개의미생물에대한설명이부족하기때문에 cell-cycle synchronization 현상을설명하기에부적절하다 (Mhaskar and Hjortso, 2002). Jones and Kompala (1999) 은유산소배양에서의 S. cerevisiae ( 또는효모 ) 의성장에관하여이미알려진대사경로를이용하여세포의진동현상을모사하였다. 이들은 glucose 발효반응식, 에탄올산화식그리고 glucose 산화식에대한반응속도식을이용하여, 세포질량, 글루코즈 / 에탄올 / 산소 / 효소농도의시간에대한변화를모델링하였다. 상기의 3 가지반응속도식에서나오는성장속도식은 batch 실험을통하여구하였으며, 여러다른모델인자들은실험결과와비교를통하여예측되었다. 이모델을통하여세포질량 / 효소농도 / 산소농도등의시간에따른진동현상을예측하였으며, 모사결과와실험결과는비교적정확하게일치하였다. Unstructured/segregated modeling ( 비체계적입자화된모델링 ) 이모델은여러상태의세포입자분포를모델링하면서비교적간단하게세포의진동현상을설명하며, 생물반응기의제어 (e.g., model predictive control, Zhu et al., 2000) 에이용하기위하여개발되었다. 이비체계적인 (unstructured) V-1

모델은세포내여러물질들의복잡한상호작용이세포성장에미치는영향을매우간략하게기술하며, 세포질량에따른여러세포들의 (segregated) 특성치를고려하기위해 population balance equation (PBE) 를이용한다. 세포입자분포식 (PBE, population balance equation) 은 rate limiting substrate 의질량수지식과결합하여사용하며, 세포내부 / 외부의복잡한대사반응식은포함하지않는다. 세포입자수지식 (PBE) 은세포성장에관한대류항 GN (convection term, ), 세포분열에의한새로운세포의생성 (division birth term, 2 pγ Ndm ), 세포분열에의한 m 어미세포의소멸 (division death term, ΓN ) 그리고 dilution 에의한세포의유실 ( DN ) 에대하여다음과같이표현된다. N( m,t ) ( G( S ) N( m,t )) m = + 2 p( m,m ) Γ ( m,s )N( m,t )dm [ D + Γ ( m )] N( m,t ) 0 t m (5-1) 여기에서, N ( m,t ) 는질량 (m) 과시간 (t) 에따른세포수, G ( S ) 는기질농도 (S ) 에따른세포성장속도, p ( m,m ) 는새롭게태어난딸세포질량 (m ) 이엄마세포의질량 m 로부터분리되어질확률함수, D 는 dilution rate, 그리고 Γ ( m,s ) 는분열강도함수이다. 분열강도함수는세포가어떤임계질량에도달함에따라세포분열을하려고하는가능성을모델링하며, 다음과같이표현된다. 0, m mt 2 ε ( m md ) Γ ( m,s ) = γ e, mt < m < md (5-2) γ, m md 여기에서 m t 는전이질량이고 m 0 는세포분열이가능하기위해서필요한어미세포의최소질량이며, m d 는분열강도값이최대인 γ 에도달할때의세포질량이다. 세포분열후새롭게태어날세포의확률함수 (newborn cell probability function) p ( m,m ) 는세포분열로부터생기는딸세포의질량분포를나타낸다. 2 2 β ( m mt ) β ( m m + mt ) α e + α e, m > m and m > mt p( m,m ) = (5-3) 0, m m and m mt 세포입자수지식 (5-1), (5-2) 그리고 (5-3) 은기질질량수지식과세포성장속도식을결합함으로서완결된다. ds = S ) D( S f S ) G( N( m,t ) dm (5-4) 0 dt Y µ ms G( S ) = (5-5) K m + S 여기에서 S 는실제기질농도, S 은 effective substrate concentration, S f 는유입되는기질농도, Y 는 constant yield coefficient, 그리고 µ m 은최대세포성장속도를뜻한다. 결론적으로, 이모델은한개의편미방식 (PDE, partial differential equation, (5-1)), 한개의미방식 (ODE, ordinary differential equation, (5-4)) 그리고한개의대수식 (AE, algebraic equation, (5-5)) 으로구성된 PDAE (partial differential algebraic equation) system 으로볼수있으며, 시간과질량을독립변수로하고 N ( m,t ), S ( t ) and G ( S ) 를종속변수로하고있다. 이식을풀기위해서는 N( m,t ) 과 S ( t ) 에대한초기조건과 N ( m,t ) 에대한경계조건이요구되며, 일반적으로이식의해석해는존재하지않고, 편미방의수치해석방법을통한수치해를얻어야한다. 이모델의모델인자는 Zhu et al. (2000) 에서찾아볼수있으며, 그림 5-1 은이모델의시간과세포질량에대한세포수의진동현상을보여준다. 그림 5-1. Numerical solution of the simplified cell PBE model (Zhu et al., 2000). V-2

Systems biology and engineering (SBE) 임영일, 년 12월 월 29일 일 임영일 2003년 이 PBE 모델에 대한 가장 간단한 수치해석법은 식 (5-1) 의 오른쪽 항들을 모두 질량방향으로 구간화 (discretization) 하여 시간에 따른 미분방정식의 집합체로 만들고 이렇게 얻어진 미방을 시간적분기 (time integrator) 를 이용하여 푸는 것이다. 이방법은 method of lines (MOL) 이라고 알려져 있고, 구간화방법에 따라서 FDM (finite difference method), FEM (finite element method) 그리고 FVM (finite volume method) 등으로 구분되어질 수 있다 (Lim et al., 2001; Mantzaris et al., 2001). Motz et al. (2002) 는 최근에 개발된 편미방 해석법인 CE/SE (conservation element and solution element) method 를 사용하여 이러한 입자수지식에 대한 매우 정확하고 빠른 수치해를 얻었다. 상기의 간단한 PBE 모델은 세포외부 환경에 대한 좀더 체계적인 모델링으로 보완된다. Mhaskar and Hjortso (2002) 는 Yeast 의 glucose-limited growth environment 하에서 세포외부 기질 농도 (glucose, ethanol), 액상/기상 산소농도, 액상/기상 이산화탄소 농도식, 그리고 세포외부에서 내부로의 글루코즈와 에탄올의 섭취 속도식을 추가하여, 즉 8 개의 미분방정식을 식 (5-1) 에 추가하여 cell PBE 를 모사하였다. 이들 추가된 8 개의 미분방정식들은 글루코즈 발효 속도식, 글루코즈 산화 속도식 그리고 에탄올 산화 속도식을 내포하고 있으므로, 세포성장속도는 세포내 대사과정에 기인하고 있다. 즉, 세포성장속도식은 세포내 글루코즈/에탄올 농도 ( G, E ) 와 세포외 액상 산소농도 ( O ) 에 따른 3 개 대사과정 (글루코즈 발효, 글루코즈 산화 그리고 에탄올 산화) 속도식의 합으로 표현된다. µ mgf G µ mgo G µ meo E K inhibit O O (5-6) + + G( S ) = K mgf + G K mgo + G K mgd + O K meo + E K med + O K inhibit + G 따라서 본모델은 1 개의 PBE 8 개의 미방 그리고 1 개의 대수식으로 기술된다. 모델인자값은 Mhaskar and Hjortso (2002) 에서 찾아볼 수 있으며, dilution rate (D=0.18 h-1) 에서 산소, 이산화탄소, 글루코즈, 에탄올 농도들 그리고 세포수 농도에 대한 진동현상을 그림 52 에서 볼 수 있다. 이 그림에서는 실험값과 모사값이 비교되어지고 있고, 비교적 서로 잘 일치하고 있음을 볼 수 있다. 이모델은 측정가능한 세포외부 물질들의 농도변화식을 추가함으로서 실험적으로 그리고 정량적으로 모델식을 확인할 수 있다. 또한 모델 인자값 추정을 위해서 실험결과를 직접 사용할 수 있는 장점을 갖는다. 그림 5-2. Comparison of plant and model concentration oscillations in continuous culture of S. cerevisiae (Mhaskar and Hjortso, 2002). Structured/segregated modeling (체계적 체계적 입자화된 모델링) 모델링 이 모델은 single-cell model 을 세포내/외부의 물질에 관하여 PBE 로 확장시킨 것이다. 즉, 시간과 내부 상태변수들 (또는 metabolites, xi, i = 1,..., L, 여기에서 L 은 세포내부 물질의 수) 을 독립변수로 하고, 이들 물질들의 질량비율 ( ψ ( x,t ) ) 을 종속변수로 하고있으므로, (L+1) 차원의 1 개 PBE 로 표현된다 (Nielsen and Villadsen, 1992). L ( G ( x, t ) ψ ( x, t )) ψ ( x,t ) i = + R( x,t ) R ( t ) ψ ( x,t ) t xi i =1 [ 여기에서, Gi ( x,t ) 는 species-i 의 생성속도이고, ] R( x,t ) 는 세포성장속도을 뜻하며, R( t ) (5-7) 는 평균 세포성장속도이다. L P R( x, t ) = γ i, j r j ( x, t ) (5-8) i =1 j =1 여기에서, P 는 세포내 대사반응식의 수이고, γ i, j 는 jth-reaction 에서의 ith-species 의 화학양론값 (stoichiometric value) 이며, r j ( x,t ) 는 대사반응속도식이다. 평균 세포성장속도 R ( t ) 는 R( x,t ) 을 모든 물질들에 대하여 적분하여 얻을 수 있다. R ( x,t ) = 0 R( x,t )ψ ( x,t )dx (5-9) V-3

편미방으로표현된식 (5-7) 에서모든내부상태변수 L 개를 K 수만큼구간화한후생기는미분방정식수 ( N ODE ) 는다음과같이표현된다. L N ODE = K (5-10) 만일 6개대사물질 (L=6) 에대하여, 6 개의대사반응식 (P=6) 을사용하고, 이들 6개의상태변수 ( 대사물질농도 ) 6 를각각 10 개씩의구간 (K=10) 으로구간화하면식 (5-7) 은 10 개의미분방정식으로정리된다. 일반적으로 15 대사물질의수가매우많은데, 15 개정도로만사용하더라도 10 개의미방이생기므로이모델링은현실적으로계산이불가능해진다. 이러한단점을보완하기위해다음과같은 cell ensemble model 이하나의대안으로여겨질수있다. 5.2 세포앙상블모델 (cell ensemble model) 세포입자수지식 (cell PBE) 은앞서지적하였듯이단일세포의성장과분열에관련된생화학적대사반응에대한설명이부족하다. 또한여러상태변수에대한 PBE 모델 (e.g., structured/segregated models) 은풀어야될너무많은미분식으로인하여매우강력한계산기를필요로한다. 이러한단점을보완하고, 단일세포내모든대사 / 조절 / 신호전달과정에대한 single-cell model (SCM) 을이용하여, 이단일세포들의집합체의동적거동을예측하기위해다음세단계로 cell ensemble model 을구성한다. i) 단일세포내여러 metabolite 들의대사과정에대한시간에따른농도변화를표현하는비선형미분방정식들 (SCM, single-cell model) 을구성한다. ii) 비균일세포군집들의동적거동을표현하기위해, 기준값에서일정하지않게분포된초기값이나대사반응상수값을갖는충분히큰수 ( 약 1000 5000) 의미생물세포들로구성된한집합체 (ensemble) 를형성시킨다. 여기에서, 비선형미분방정식의초기값과대사반응상수값들의 random perturbation 은평균값에서실험적으로알려진일정한분산값을갖는정규분포를이룬다고가정한다. iii) 세포외부의먹이나배설물질에대한농도변화를표현하는미분방정식을추가한다. 따라서이모델식은다음과같이표현된다. dxi = fi( xi, θi,y ), 1 i M dt (5-11) dy = g( x1,...,xn,y, φ ) dt 이식은세포내여러물질들의농도를나타내는상태변수 (x), 세포내대사반응상수 (θ), 세포외부물질에관한상태변수 (y), 세포외부대사반응에관련된반응상수및전달계수 (φ ) 그리고세포수 (M) 등으로표현된다. 이모델은 Ataai and Shuler (1985) 에의하여처음으로개발되었고, 최근에 Henson et al. (2002) 는효모 (yeast) 에대하여 1000 개의세포로구성된집합체를사용, glycolytic oscillation 을설명하였다. 그림 5-1 은 yeast 가세포외부로부터 glucose 를섭취하고, 여러효소들의촉매작용으로 glycerol 과 ethanol 을형성하며, 부산물을배출하는무산소글루코즈대사 (anaerobic glycolysis metabolism) 반응경로를보여준다. 그림 5-3. Yeast glycolytic reaction pathway for anaerobic growth (Henson et al., 2002) V-4

세포벽을통과한 glucose (S1) 가두개의 ATP 를소비하여 2 개의 G3P/DHP (S2) 를만들고 (v1), NADH 1 개를소비하여 glycerol 를만들며 (v6), NADH 1 개를생성하면서 1,3-BPG (S3) 로변환된다 (v2). 반응 v3 는 2 개의 ATP 를생성시키면서 1,3-BPG 로부터 acetaldehyde/pyruvate (S4) 로변환시키고, 이 S4 는 1 개의 NADH 를소비하면서 ethanol 을만든다 (v4). acetaldehyde/pyruvate (S4) 는또한세포벽을통과하여외부로배출되어또다른물질로변환된다 (v7). v3 와 v5 는 energy metabolism 의반응속도식을단순하게표현한다. 그림 5-3 로부터얻은대사경로와반응속도식 (v1-v6) 으로부터 6 개의미분방정식을그림 5-4 에서볼수있다. 이식에서 intermediate metabolites (S1-S4) 와 ATP (A 3 ), NADH (N 2 ) 에대한시간에따른농도변화는알려진반응속도상수값 (k 1 -k 6 ) 과 glucose 섭취속도값 (J 0 ), 그리고주어진여러다른모델인자값을이용하여계산할수있다 (Henson et al., 2002). 그림 5-4. Yeast glycolytic intracellular reactions. 이들 6 개의식에세포외부로배출되는 acetaldehyde/pyruvate (S4 ex ) 에관한한개의미분방정식을추가하여간략화된 Yeast glycolytic reaction 을완결한다. ex ds M M 4 ϕ ϕ ex = Ji v7 = Ji ks4 (5-12) dt M 1 M 1 여기에서 M 은세포수를뜻하며, Henson et al. (2002) 는 1000 개의세포에관하여 acetaldehyde/pyruvate (S4) 의 ex 배출속도값 Ji = κ ( S4,i S4 ) 을이용하여총 6001 (= 6 1000 + 1 ) 개의미분방정식을수치적으로해석하였다. 비균일세포상태모사를위해, i) 세포내중간체들 (S 1,i -S 4,i for i=1 M) 과 ATP, NADH 농도 (A 3,i, N 2,i for i=1 M) 의초기값 ( t = 0 ) 이평균값에서일정한분산값을갖는정규분포를이룬다고가정하거나, ii) 세포내대사반응상수값 (k 1,i -k 6,i for i=1 M) 이평균값에서일정한분산값을갖는정규분포를이룬다고가정한다. 그림 5-5 는 6000 개미분방정식의초기값을 random perturbation 시킨후얻은모사결과를보여준다. 초기에비동기화된세포집합체들은시간이지남에따라, 세포외부의 acetaldehyde/pyruvate (S4 ex ) 농도의중첩효과 (coupling effect, see 식 (5-12)) 로인하여, NADH 농도가점차적으로동기화 (synchronization) 되면서 25 분후부터는 1 분간격으로진동현상을보여주고있다. 이모델은 structured/segregated modeling 기법과맥을같이하고있으며, structured/unsegregated modeling 기법에서확장된것으로볼수있다. cell ensemble model 은본보고서의 3장에서기술한세포모사에관한모델식을직접사용할수있지만, 대사반응식의수가많아지면서풀어야할미분방정식수 ( H M, 여기서 H 은대사반응식수그리고 M 은세포수 ) 가증가하므로, 계산시간을단축하기위해서는적절한수의선택이요구된다. 그림 5-5. Cell population synchronization for a 1000-cell ensemble where the initial conditions are perturbed by a Gaussian distribution (Henson et al., 2002). V-5

5.3 결론 Flow cytometry 는각각의세포내물질들 (proteins, DNA etc.) 의분포를측정할수있는매우강력한방법으로부상하고있다. 또한 Flow cytometry 를이용하여최근에온라인으로세포의특성치에대한분포를측정할수있게되었다. 이것으로부터얻은실험결과를바탕으로세포군집들의동적거동에대한모델링은생물공정의설계 / 운전 / 제어에큰도움을줄것이다. 비균일세포군집들의동적거동과분포에관한모델링은크게입자수지식 (PBE, population balance equation) 과세포앙상블모델 (cell ensemble model) 로나누어볼수있다. 입자수지식은세포내대사반응을매우간략화시킨후세포질량에따른세포수의분포를기술하는방법 (unstructured/segregated modeling) 과세포내대사반응을포함하면서이들의 population balance 를고려하는방법 (structured/segregated modeling) 등이있는데, 현실적인계산성능을고려해볼때, unstructured/segregated modeling 이많이사용되고있다. 세포앙상블모델은복잡한대사반응을갖는단일세포들의집합체로서 cell population dynamics 를설명하며, 단일세포모델 (SCM, single-cell model) 의 random perturbation 을바탕으로생물반응기내수백만의세포군집들의분포를예측한다. PBE 모델은 cell-cycle dependent oscillation 을예측하기위하여, 그리고, cell ensemble model 은단일세포내 glycolytic oscillation 을설명하는데주로이용되었으며, 비록두모델은다른접근방식을취하고있지만, 비균일세포들의동적거동을예측할수있다. References Ataai, M. M. and M. L. Shuler (1985), Simulation of CFSTR through development of a mathematical model for anaerobic growth of E. coli cell population, Biotech. Bioeng., 27, 1051-1055. Henson, M. A. (2003), Dynamic modeling of microbial cell populations, Current Opinion in Biotechnology, 14, 460-467. Henson, M. A., D. Muller and M. Reuss (2002), Cell population modeling of yeast glycolytic oscillations, Biochem. J., 368, 433-446. Jones, K. D. and D. S. Kompala (1999), Cybernetic model of the growth dynamics of S. cerevisiae in batch and continuous cultures, J. Biotechnol., 71, 105-131. Lim, Y. I., J. M. Le Lann & X. Joulia (2001), Accuracy, temporal performance and stability comparisons of discretization methods for the solution of Partial Differential Equations (PDEs) in the presence of steep moving fronts, Comp. Chem. Eng., 25(11-12), 1483-1492. Mantzaris et al. (2001), Numerical solution of multivariable cell population balance models: I, II and III, Computers and chemical engineering, 25, 1411-1481. Mhaskar, P., M. A. Hjortsø and M. A. Henson (2002), Cell population modeling and parameter estimation for continuous cultures of S. cerevisiae, Biotechnol. Prog., 18, 1010-1026. Motz, S., A. Mitrovic & E.-D. Gilles (2002), Comparison of numerical methods for the simulation of dispersed phase systems, Chem. Eng. Sci., 57, 4329-4344. Nielsen, J. and J. Villadsen (1992), Modeling of microbial kinetics, Chem. Eng. Sci., 47, 4225-4270. Zhang, Y., A. M. Zamamiri, M. A. Henson and M. A. Hjortsø (2002), Cell population models for bifurcation analysis and nonlinear control of continuous yeast bioreactors, J. Process Control, 12, 721-734. Zhu, G.-Y., A. Zamamiri, M. A. Henson and M. A. Hjortsø (2000), Model predictive control of continuous yeast bioreactors using cell population balance models, Chem. Eng. Sci., 55, 6155-6167. V-6