1. sports
2. skating
3. figure skating

Mathematical Modeling of Aquaporin-4 Intracellular Transport and Complex Column
Distillation Process
SEON B. KIM
B.S., Yonsei University, Seoul, Korea, 2004
THESIS
Submitted as partial fulfillment of the requirements for the degree
of Doctor of Philosophy in Bioengineering
in the Graduate College of the
University of Illinois at Chicago, 2013
Chicago, Illinois
i
This dissertation is dedicated to my parents and wife, Taeyeon, without whom it would never
have been accomplished.
ii
ACKNOWLEDGEMENTS
I would like to express my deepest gratitude to my advisor, Prof. Andreas Linninger, for
his excellent guidance, support, and providing me with an invaluable atmosphere for doing
research. I am also grateful to my committee members Prof. Christos Takoudis, Prof. Umila
Diwekar, Prof. Seungpyo Hong, and Prof Simon Alford for their insights and suggestions which
added great value to this work.
I would also like to thank my former graduate student colleagues at LPPD who have
helped in various aspects of this work: Dr. Gerald Ruiz for the discussion and advice for the
distillation project, Dr. Madhawa Hettiarchchige for his special assistance in fluid transport
visualization, Jeonghwa Moon in the optimization problem and programming technique, and
Ying Hsu for the sharing of her knowledge in the AQP4 transport project. I would also like to
acknowledge the valuable collaboration with LPPD graduate students; Dr. Sukhraaj Basati, Brian
Sweetman, Eric Lueshen, Chih-Yang Hsu, Kevin Tangen, Thomas Marrinan, Indu Venugopal, as
well as undergraduate researchers Minh Tran, Cierra Hall, Nairyna Constantino, Maurice
Chojechi, Masood Qader, Tejen Soni, Joe Lancaster, Bhargav Desai, Fabian Davalos, Sebastian
Pernal, Sabarish Chocklanigean.
I would also like to thank my parents for their devoted support, and elder sister and elder
brother. They were always supporting me and encouraging me with their best wishes.
I would like to thank to my son and daughter. Their delightful welcome at home was
always very precious to me and cheering me up to finish this work.
iii
TABLE OF CONTENTS
CHAPTER
PAGE
Part I. Stochastic Reaction-Diffusion Model of AQP4 Intracellular Polarization............................ 1
1
LITERATURE REVIEW OF AQUAPORIN-4 SYNTHESIS AND INTRACELLULAR
TRANSPORT ................................................................................................................................... 3
1.1
Discovery of Extracellular Water Transport Channels ................................................... 3
1.2
Aquaporin-4 (AQP4) in the Central Nervous System ..................................................... 4
1.2.1
Cytotoxic edema ........................................................................................................ 4
1.2.2
Vasogenic edema ....................................................................................................... 5
1.3
A Morphology of Astrocyte Intracellular Space ............................................................. 5
1.3.1
The Pathological Role of Astrocytes ........................................................................ 6
1.3.2
Water Homeostasis Functionality of Astrocyte ....................................................... 7
1.3.3
The Organization of Microtubules ............................................................................ 8
1.4
AQP4 Synthesis Mechanisms ........................................................................................... 9
1.4.1
Transcriptional Activation by Transcription Factor, Nrf2 .................................... 10
1.4.2
AQP4 Transcription ................................................................................................. 11
1.4.3
AQP4 Translation .................................................................................................... 11
1.4.4
Hydrocephalus Inducing AQP4 Up-regulation ...................................................... 12
1.4.5
Sulforaphane Enhancing Aquaporin-4 Expression ................................................ 13
iv
1.5
2
Aquaporin-4 Transport Mechanisms .............................................................................. 14
METHODOLOGIES
ON
STOCHASTIC
SIMULATION
OF
AQUAPORIN-4
INTRACELLULAR TRANSPORT.............................................................................................. 17
2.1
Intracellular Structure and Conditions ........................................................................... 19
2.1.1
Aquaporin-4 Water Channels in Cell Membrane .................................................. 21
2.1.2
Modeling of Microtubule Structure ........................................................................ 23
2.2
Stochastic Model and Simulation of AQP4 Synthesis .................................................. 27
2.2.1
Kinetic Reaction Model of AQP4 Synthesis .......................................................... 27
2.2.2
Transcription Kinetics ............................................................................................. 31
2.2.3
Translation Kinetics ................................................................................................. 32
2.2.4
Degradation Kinetics ............................................................................................... 32
2.3
Stochastic Simulation of the Intracellular Transports ................................................... 33
2.3.1
Mechanisms of AQP4 Transport............................................................................. 34
2.3.2
The Chemical Master Equation ............................................................................... 37
2.3.3
Monte Carlo simulation and the Gillespie algorithm ............................................ 39
2.3.4
Gillespie Multi-Particle (GMP) Method with Modified Diffusion Fraction for the
Structured Coordinate System ................................................................................................ 40
2.3.5
Gillespie Multi-Particle (GMP) Method Event Time Computation for the
Unstructured Coordinate System ........................................................................................... 42
v
3
RESULT OF THE AQP4 INTRACELLULAR TRANSPORT BY STOCHASTIC
SPATIOTEMPORAL REACTION-DIFFUSION SIMULATION ............................................ 44
3.1
Astrocyte Cell Reconstruction ........................................................................................ 44
3.2
Microtubule Compartment Implementation .................................................................. 45
3.3
Kinetic Model of AQP4 Synthesis ................................................................................. 47
3.3.1
AQP4 Synthesis at Steady State.............................................................................. 47
3.3.2
Dynamic AQP4 Synthesis ....................................................................................... 48
3.4
Spatial Translocation of Aquaporin-4 ............................................................................ 50
3.5
Nrf2 activation by SFN treatment in astrocytes ............................................................ 52
3.6
Up-regulation of AQP4 by SFN .................................................................................... 53
4
DISCUSSION.......................................................................................................................... 55
5
CONCLUSION ....................................................................................................................... 58
CITED LITERATURE....................................................................................................................... 59
Part II. Design and Modeling of Energy Efficient Complex Distillation Process .......................... 65
1
INTRODUCTION
TO
THE
DESIGN
AND
SYNTHESIS
OF
COMPLEX
DISTILLATION PROCESS.......................................................................................................... 67
2
METHODOLOGIES OF THE COMPUTATIONAL MODELING OF COMPLEX
COLUMN CONFIGURATION .................................................................................................... 70
2.1
Vapor-Liquid Phase Equilibrium Model ........................................................................ 70
vi
2.1.1
Constant Relative Volatility Model ........................................................................ 70
2.1.2
Raoult’s Law Model ................................................................................................ 71
2.2
Complex Column Model ................................................................................................. 72
2.2.1
Component Balances for Complex Columns ......................................................... 73
2.2.2
Continuous Model for Complex Columns ............................................................. 78
2.3
Temperature Collocated Continuous Profile Equation ................................................. 79
2.4
Computation of Composition Profiles............................................................................ 80
2.4.1
2.5
Robust Pinch Point Computation ................................................................................... 83
2.5.1
Pinch point location for constant relative volatilities ............................................ 84
2.5.2
Pinch Point Location for for Ideal Solutions .......................................................... 86
2.5.3
One-dimensional deflation method for Pinch Point Calculation .......................... 88
2.5.4
Pinch Point Location for for Nonideal Solutions ................................................... 89
2.6
Piecewise polynomial approximation ............................................................................ 90
2.7
Inverse Design Problem .................................................................................................. 92
2.7.1
3
Computation of Nonideal Mixture Model .............................................................. 82
Bubble Point Distance (BPD) for Feasible Design ................................................ 93
RESULTS OF DESIGN AND SYNTHESIS OF COMPLEX COLUMN DISTILLATION
NETWORKS .................................................................................................................................. 97
3.1
Ternary Mixture Separation with Complex Column Networks ................................... 97
vii
3.2
Quaternary Ideal Mixture Separation with Complex Column Networks ...................102
3.2.1
Case study I – Equimolar Alkane Mixture Separation ........................................103
3.2.2
Case study II – Non-equimolar Alkane Mixture Separation ...............................114
3.2.3
Case study III – Aromatic and Alkane Mixture Separation ................................117
4
DISCUSSION........................................................................................................................126
5
CONCLUSION .....................................................................................................................129
CITED LITERATURE.....................................................................................................................131
viii
LIST OF TABLES
TABLE
PAGE
Table I. Fixed compartments design specification ........................................................................... 21
Table II. Kinetic rates of AQP4 synthesis ......................................................................................... 48
Table III. Definitions of vapor flow, liquid flow, vapor and liquid compositions for all four
sections of a complex column ................................................................................................... 78
Table IV. Pseudo-code implementation to compute the liquid composition profiles using
orthogonal collocation on finite elements for each column..................................................... 81
Table V. Design and problem specifications and internal vapor and liquid flows. ........................ 98
Table VI. Design and problem specifications and internal vapor and liquid flows based on
Underwood method for the complex column configuration of Figure 30. ...........................100
Table VII. Antoine parameters of the alkane mixture of Pentane (A), Hexane (B), Heptane (C),
and Octane (D) used for the calculation of the partial vapor pressure. .................................103
Table VIII. Vapor rate for each column and total vapor rate of each designed network. ............107
Table IX. Column design aspects obtained by solving the inverse design problem of nine column
configurations to separate a quaternary mixture using the minimum BPD design algorithm.
In this table, continuous length were converted to actual tray using eq. (71); remark were
rounded to next integer for better comparison with Aspen simulation. ................................107
ix
Table X. Detailed network design specifications and operating conditions obtained by solving
the inverse design problem to separate a quaternary mixture of Pentane, Hexane, Heptane,
and Octane using two different molar compositions of equal molar FEED 1 of [0.25; 0.25;
0.25; 0.25] and FEED2 of [0.1780; 0.4046; 0.2771; 0.1403] in three different complex
networks NW1, NW2, and NW3 as in Figure 34. ....................................................................115
Table XI. Detailed network design specifications and operating conditions obtained by solving
the inverse design problem to separate a quaternary mixture of Benzene, Toluene, Octane,
and Nonane using two different molar compositions of equal molar FEED 1 of [0.25; 0.25;
0.25; 0.25] and FEED2 of [0.1780; 0.4046; 0.2771; 0.1403] in three different complex
networks NW1, NW2, and NW3 as in Figure 34. ....................................................................119
x
LIST OF FIGURES
FIGURE
PAGE
Figure 1. Organization of microtubules in various astrocytes; (a) immunofluorescence image of
microtubules of astrocytes using an anti-tunulin antibody (Tomás et al. 2005), (b)
Fluorescence speckle microscopy illuminating the intricate network f microtubules (yellow)
and actin filament (Purple) (image credit: Torsten Wittmann, UCSF). .................................... 9
Figure 2. The organization of AQP4 gene and two isoforms of aquaporin-4 (Umenishi and
Verkman 1998). .......................................................................................................................... 10
Figure 3. Schematic representations of active bidirectional transport by motor proteins along
microtubules and active filaments in a cell. Myosin motor proteins walk along actin
filaments (red) at the cortex. Microtubule-based motors including the kinesin and dynein
motor proteins walk to the plus ends of microtubules (green), which are oriented toward the
actin cortex and to the minus end of the microtubule, which is located at the microtubuleorganizing center (MTOC, green) near the cell nucleus (blue) respectively. Figure credit:
(Ross et al. 2008) ........................................................................................................................ 16
Figure 4. Sketch of intracellular structure of astrocyte. ................................................................... 19
Figure 5. Schematic of OAP assembly from AQP4 M1/M23 heterotetramers. (A) (top) M23 and
M1 AQP4 monomers. Projection from M23 monomer denotes potential site of M23-M23
intertetrameric
interaction.
(bottom)
Six
possible
independent
AQP4
M1/M23
heterotetramers. Four-way rotational variants of each are possible as well. (B) Examples of
OAPs at low and high M23:M1 ratios. (Jin et al. 2011) .......................................................... 22
xi
Figure 6. An example of Smooth and Shortest Path Finding algorithm in 2d. All paths avoid the
obstacles (rectangular shape) on the straight way to the target points (marked as ‘x’). ........ 25
Figure 7. Schematic mechanism of AQP4 protein synthesis. In this study, Nrf2 activation by
SFN treatment up-regulates AQP4 expression on the astrocyte cell. ..................................... 28
Figure 8. Cell molecular simulator flowchart. .................................................................................. 29
Figure 9. Simplified transporting cascade of protein molecule between compartments. .............. 36
Figure 10. Structured meshes in the Cartesian coordinate system for different directional
diffusion ...................................................................................................................................... 41
Figure 11. Diffusion event simulation by Gillespie multi-particle method (GMP). The fluctuation
of diffusing molecules is caused by randomness of event and discreteness of number of
molecule. In the original GMP method does not account for the empty site situation after
single diffusion event simulation when a group of molecules is initially located only in the
specific sub-volumes or where a site has all empty neighboring sites.................................... 41
Figure 12. Diffusion event simulation by modified fractional Gillespie multi-particle method.
Fractional diffusion prevents the undesired empty site problem happened at initial steps of
diffusion by source point or when a site are surrounded with all empty neighboring sites.
Modifying diffusion ratio from 1 to 0.8 (typically we use 80% of clearance of diffusion
event site) shortens simulation interval from ∆τ to 0.8*∆τ. ..................................................... 42
xii
Figure 13. Stochastic Reaction-Diffusion Flowsheet. The fractional diffusion is not considered in
this conceptual flowsheet. .......................................................................................................... 43
Figure 14. Cell reconstruction. The series of confocal microscopic images (a) are reconstructed
to the computational model (b). The computational model has 11,390 tetrahedron volume
meshes. ........................................................................................................................................ 44
Figure 15. Cell meshing post-processing of wrapping and smoothing. Selecting proper
parameters in Mimics mesh generator, we can obtain the high quality 3d model of astrocyte
cell. .............................................................................................................................................. 45
Figure 16. Simple implementation of straight microtubule structures. Microtubules share the
mesh information with cell meshes. (a) shows the test of active transport by microtubules
and (b) is the 1st cell model with consideration of microtubules volume ratio. Microtubules
facilitate the vesicle transport to the cell boundary. ................................................................. 46
Figure 17. Real-like microtubule implementation in the astrocyte 3d cell model. The red curves
represent microtubules. .............................................................................................................. 47
Figure 18. Stochastic simulation results of AQP4 synthesis. Solid lines in red color represent the
mean of 100 stochastic simulations and the black dotted lines are the result of the
deterministic model. The proposed model and kinetic parameters match the experimental
results. ......................................................................................................................................... 49
xiii
Figure 19. Stochastic active and passive transport result with different number of microtubules
during 100 simulation time in the circular 2d cell model. The more microtubules in the cell,
the more it facilitates the active transport to the boundary of cell model which is assumed
plasma membrane in the real 3d cell model. In the simulation of 1000 molecules, there are
21 molecules for the case of the number of 5 MTs, 93 molecules for 100 MTs and ~390
molecules for 300 MTs associated with active transport by MTs........................................... 50
Figure 20. AQP4 in primary cultured astrocytes (a) and model (b). In (a), AQP4 vesicles are
stained as green dot and blue circle in the center is nucleus. AQP4 are seen in vesicles as
well as along astrocytic processes during their transport towards the endfeet. Scale bar in
both images indicates 10 microns. The computational model of astrocyte in (b) visualizes
realistic microtubules and size-varied blue dot AQP4 vesicles in 3d distributed system...... 52
Figure 21. Translocation of Nrf2 from cytoplasm to nucleus after SFN exposure. In control cells,
Nrf2 is found in cytoplasm with a weak basal immunoreactivity in the nucleus. After SFN
exposure, Nrf2 is translocated to the nucleus showing a high immunoreactivity, and Nrf2filled vesicles marking transport of Nrf2 are observed. ........................................................... 53
Figure 22. Western blotting shows that SFN induced a 48% up-regulation after 9 hours of
continuous exposure, and 68% up-regulation after 18 hours compared to control. ............... 54
Figure 23. A complex column network with two feed streams and three product streams. (a)
depicts skeptical configuration of simple and complex column configuration equivalent to
Petlyuk column and (b) shows the sectional decomposition of complex column with
xiv
implemented nomenclature of vapor (Vi), liquid (Li), concentration (xi), and Feed (Fj) for
each section i or feed input j. ..................................................................................................... 72
Figure 24. Mass Balance at the top section S1 in Figure 23 (b). ..................................................... 75
Figure 25. Mass balance at the bottom section S4 in Figure 23 (b). ............................................... 76
Figure 26. A column profile map for a mixture with constant relative volatilities of (A=5, B=2,
C=1) with only one real solution at a reflux of 0.6 and XΔ= [1.2, -0.1, -0.1]. The thick black
trajectory indicates a realizable liquid composition profile of a typical column section....... 86
Figure 27. A one dimensional pinch point temperature search showing all pinch temperatures
where f(T)=0 for XΔ=[0.9, 0.05, 0.03, 0.02] and RΔ=5, accompanied by a table showing all
pinch temperatures and compositions. The pinch temperatures indicated on the graph (T1T4) are the bubble point temperatures of the four pinch points. .............................................. 87
Figure 28. (a) Surrogate function g1(T) vs. bubble point temperature where the highest pinch
point has been eliminated, (b) surrogate function g2(T) vs. bubble point temperature where
the two highest pinch points have been eliminated, and (c) surrogate function g2(t) vs.
bubble point temperature where the three highest pinch points have been eliminated. ........ 89
Figure 29. Detailed information flow diagram to find feasible design specifications for an entire
flowsheet. Stage 1 uses the feed and sets up the final product purities. In stage 2, initial
guess of intermediate products purities and operating conditions (reflux or reboil ratios) are
specified. Then product flowrate, internal vapor and liquid flowrate, stationary points and
xv
difference points parameters of each column section are computed using species balances.
Stage 3 performs BPD criterion with obtained composition profiles by OCFE and
polynomial arithmetic to assess feasibility. The schematics on the right side depict the
specifications for the computations of each stage. ................................................................... 96
Figure 30. Complex column arrangement to separate a ternary mixture ........................................ 97
Figure 31. Column profiles for a complex column network shown in Figure 30. ......................... 99
Figure 32. Complex column arrangement to separate a ternary mixture ......................................101
Figure 33. 18 Separation flowsheets to separate a quaternary mixture. These configurations
contain simple columns or basic complex columns. ..............................................................103
Figure 34. 9 optimal basic network configurations to separate a quaternary mixture. These
configurations have an optimal tradeoff between capital and operating costs, and minimum
energy requirement...................................................................................................................104
Figure 35. Frame (a) shows the flowsheet with feeds (F1-F4) and products (P1-P4) for the
separation of a quaternary mixture. Frames (b) to (d) represent the tetrahedral composition
profiles of columns C1–C3 for the four species A, B, C, and D. The blue trajectories stand
for rectifying sections and red trajectories for stripping section. ..........................................105
Figure 36. Initialization of AspenPlus by temperature collocation approach to validate network
NW1. Frames (a) and (b) show the composition profiles of all columns in the stage number
and temperature space respectively found by temperature collocation. Frames (c) and (d)
xvi
describe the profiles of the entire network in the stage number and temperature space
respectively found by AspenPlus. (remarks: component A: blue,--x--; component B: green,-●--; component C: pink,--■--; component D: red,--+--) ........................................................109
Figure 37. Initialization of AspenPlus by temperature collocation approach to validate network
NW2. Frames (a) and (b) show the composition profiles of all columns in the stage number
and temperature space respectively found by temperature collocation. Frames (c) and (d)
describe the profiles of the entire network in the stage number and temperature space
respectively found by AspenPlus. (remarks: component A: blue,--x--; component B: green,-●--; component C: pink,--■--; component D: red,--+--) ........................................................110
Figure 38. Initialization of AspenPlus by temperature collocation approach to validate network
NW3. Frames (a) and (b) show the composition profiles of all columns in the stage number
and temperature space respectively found by temperature collocation. Frames (c) and (d)
describe the profiles of the entire network in the stage number and temperature space
respectively found by AspenPlus. (remarks: component A: blue,--x--; component B: green,-●--; component C: pink,--■--; component D: red,--+--) ........................................................111
Figure 39. Composition profiles versus number of stages and temperature stage domains found
by temperature collocation for networks NW4, NW5, and NW6. Frame (a) shows the
composition profiles of all columns in the stage number and temperature space for NW 4,
NW5, and NW6. Frame (b) contains the schematic flowsheet diagrams for all three
networks. (remarks: component A: blue,--x--; component B: green,--●--; component C:
pink,--■--; component D: red,--+--) ........................................................................................112
xvii
Figure 40. Composition profiles versus number of stages and temperature stage domains found
by temperature collocation for networks NW7, NW8, and NW9. Frame (a) shows the
composition profiles of all columns in the stage number and temperature space for NW 7,
NW8, and NW9. Frame (b) shows the schematic flowsheet for each network. (remarks:
component A: blue,--x--; component B: green,--●--; component C: pink,--■--; component
D: red,--+--) ..............................................................................................................................113
Figure 41. Liquid composition profiles on each column of the networks NW1, NW2, and NW3.
These feasible designs were found by temperature collocation method. The quaternary
mixture MIX2 constituted by Pentane (A), Hexane (B), Heptane (C), and Octane (D) is
separated in each one of the three networks with feed, FEED 2 of [0.1780; 0.4046; 0.2771;
0.1403]. The composition profiles on the top are corresponding to first column, second row
to second column, and bottom profiles to third column for each network. In the composition
triangle, the labels, Fi, refer to the feed tray compositions, not the feed compositions. ......116
Figure 42. Total vapor flow rate required for the separation of alkane quaternary mixture of
Pentane, Hexane, Heptane, and Octane using two different molar compositions of equimola
feed of [0.25; 0.25; 0.25; 0.25] and non-equimolar feed of [0.1780; 0.4046; 0.2771; 0.1403]
in three different complex networks NW1, NW2, and NW3. .................................................117
Figure 43. Liquid composition profiles on each column of the networks NW1, NW2, and NW3.
These feasible designs were found by temperature collocation method. The quaternary
mixture of Benzene, Toluene, Octane, and Nonane is separated in each one of the three
networks with equimolar feed, FEED1, of [0.25; 0.25; 0.25; 0.25]. The composition profiles
xviii
on the top are corresponding to first column, second row to second column, and bottom
profiles to third column for each network. In the composition triangle, the labels, Fi, refer to
the feed tray compositions, not the feed compositions. .........................................................120
Figure 44. Liquid composition profiles on each column of the networks NW1, NW2, and NW3.
These feasible designs were found by temperature collocation method. The quaternary
mixture of Benzene, Toluene, Octane, and Nonane is separated in each one of the three
networks with non-equimolar feed, FEED2 of [0.1780; 0.4046; 0.2771; 0.1403]. The
composition profiles on the top are corresponding to first column, second row to second
column, and bottom profiles to third column for each network. In the composition triangle,
the labels, Fi, refer to the feed tray compositions, not the feed compositions......................121
Figure 45. Total vapor flow rate required for the separation of aromatic and alkane quaternary
mixture of Benzene, Toluene, Octane, and Nonane using two different molar compositions
of equimola feed of [0.25; 0.25; 0.25; 0.25] and non-equimolar feed of [0.1780; 0.4046;
0.2771; 0.1403] in three different complex networks NW1, NW2, and NW3. ......................122
xix
SUMMARY
Mathematical modeling is a fundamental approach for many engineering problems,
which enables to solve many complex transport phenomena and bioengineering problems thus
far classical approach could not solve. Moreover, powerful computational performance using
improved computer technologies may look up the solutions obtained by traditional methods to be
more accurate and appropriate to the aim of the problem. For example, the classical design
procedure of the distillation system is not suitable because the design variable such as total
number of tray, reflux ratio, non-key component composition, etc. must be optimized through the
rigorous computational simulation.
In this thesis, we focused on the mathematical model and its solution to find more
accurate and better procedure comparing to other classical methodologies for both
bioengineering and chemical engineering systems. Project we have investigated in
bioengineering system is AQP4 intracellular transport and we tried to obtain computational
AQP4 protein distribution to the target areas. For the chemical engineering project, we
researched distillation process with complex column to find the optimized operating condition to
save energy. Each problem used mathematical modeling to describe unprecedented visualization
of aquaporin-4 intracellular transport and the design and synthesis methodologies of the complex
column configuration to separate up to 4 component petroleum mixtures respectively.
This thesis is outlined with two project-based parts. In part I, the problem of aquaporin-4
(AQP4) intracellular transport is addressed. Integration of biological experiment and
computational result enhanced the accurate expectation of AQP4 regulation and its 3d visual
xx
outcome made a clear delivery of the messages on the AQP4 transport and regulation
mechanism. The stochastic model and simulation for the AQP synthesis and transport coped with
the uncertainty of the cellular level mechanism and its validation with deterministic model and
experimental results enhanced the validity of the proposed mathematical model. The presented
approach will be beneficial to expect the amount of up-regulated AQP4 and find the control
mechanism of AQP4 concentration and also look at the possible variation of the in-vivo test.
Part II addresses computational solution of the design of complex column distillation
process in terms of energy efficiency. Rigorous simulation with different 4 component mixture
system for the 18 complex column configurations was shown in this part. The available
operating parameters includes all operating conditions to be considered in the optimized design
procedure and conclusive composition profiles were obtained satisfying the design criteria such
as product purity. To accomplish this project, the numerical approach to solve pinch points, less
dimensional piece-wise composition profiles, and inverse design for the minimum bubble point
distance was carried out. The conclusive messages that one of complex column networks to
separate different mixture systems is best efficient configuration with obtained operating
condition as well as that a different mixture system may have a different best configuration were
derived from this study. This research will benefit the remediation of the classical simple column
distillative separation in terms of energy efficiency.
xxi
PART I. STOCHASTIC REACTION-DIFFUSION MODEL OF AQP4
INTRACELLULAR POLARIZATION
Water homeostasis in a cell is now the center of physiological importance in many brain
pathologies. Many types of brain edema are revealed to be associated with mortality and
morbidity, and brain swelling increases intracranial pressure, and magnitude of cerebral edema is
now known proportional to the quantity of water in the brain. The mechanisms underlying the
development of water imbalance is still unclear, however, recent studies indicate that aquaporins
(AQPs) in the astrocytes of which processes envelope synapses of neurons may have potential
roles in the water transport and accordingly the formation of brain edema (Badaut et al. 2002;
Saadoun et al. 2002). Among discovered 13 aquaporin families (Verkman 2011), the named
aquaporin-4 (AQP4) is most abundant protein in the brain cell and its activity is being studied as
a critical element on the water transport regulation (Frigeri et al. 1995; Venero et al. 2001;
Tomás-Camardiel et al. 2004; Finnie et al. 2008).
In this part, we aimed to the synthesis and intracellular transport of AQP4 to understand
of its regulation in responding to increased transcription factor, Nrf2, expression through the
epigenetic alteration such as post-injury sulforaphane (SFN) administration. We proposed the
mathematical modeling of the AQP4 protein synthesis and spatiotemporal distribution of upregulated AQP4 into cytoplasm targeting to the end-feet of astrocyte processes. We used one of
stochastic simulation algorithms to demonstrate that the dynamics of the model is to close the
nature of living cells. In the synthesis of AQP4, the biochemical reaction is the main part of
mathematical modeling, in which sophisticated scenarios such as independent or sequential
interactions take place. Once proteins secreted out of nucleus in the form of vesicles, it travels to
1
2
other cell compartments via two separate modules of transport; passive transport by free
diffusion through the cytoplasm and active transport by motor proteins along microtubules.
With biological understanding of transcription and translation of AQP4, we set a
mathematical model using propensity function which is identical to the reaction rate equation of
the ordinary differential equations (ODEs) in the deterministic approach and simulated it with
time discretized approximation of Gillespie Algorism which is the most popular stochastic
algorithm. Secreted AQP4 protein transport from nucleus pores through endoplasmic reticulum
(ER) to the plasma membrane (PM) is modeled by stochastic diffusion-reaction and deterministic
convection along microtubules. The spatial intracellular domain was created from mesh
generator, Matrialize MIMICS v15 (Materialise), using 40 series of confocal microscopic images
cultured in Dr. Alford’s lab at UIC and the meshed volumes were taken post-processing of zone
separation by ANSYS ICEM-CFD (Ansys 2013) for the computational purposes. Generic
structures of microtubules created by shortest and smooth path algorithm were implemented in
the spatial domain and the spatiotemporal simulation of AQP4 polarization was carried out with
engineering simulation tool, MATLAB v2012 (MathWorks 2013).
This part is organized with following structures. Chapter 1 reviews topological biology of
astrocyte and AQP4 as well as considered intracellular functionalities. Chapter 2 introduces
methods of mathematical modeling used in this study and simulation of intracellular AQP4
polarization for the proposed case studies. Results and visualization of simulation follows in
chapter 3 and then lastly, this part is closed with discussion in chapter 4 and conclusion in
chapter 5.
3
1
LITERATURE REVIEW OF AQUAPORIN-4 SYNTHESIS AND
INTRACELLULAR TRANSPORT
1.1
Discovery of Extracellular Water Transport Channels
The existence of a water channels on the cell membrane was realized prior to the
identification of the protein, based on observations that Red Blood Cells (RBCs) have higher
water permeability than would be expected from an equivalent area of bilayer lipid membrane
(Benga et al. 1986; Denker et al. 1988; Preston et al. 1992). The first water channel protein was
discovered in human RBC membranes by Benga et al. (Benga et al. 1986) and it was first
purified by Denker et al. (Denker et al. 1988). Provisionally it was called CHIP28 (Preston et al.
1992) and later named aquaporin-1 (AQP1) by Agre et al. (Agre et al. 1993).
To date, at least 13 AQPs, numbered 0 to 12, have been identified in various mammalian
tissues (Verkman 2011). AQP3, AQP7, AQP9 and AQP10 have been shown to transport nonionic small solutes, such as urea and glycerol, in addition to water, whereas AQP1, AQP2, AQP4,
AQP5 and AQP8 are highly selective to water and exclude small solutes. AQP6 has been
identified as an intracellular vesicle water channel with anion permeability that is activated by
low pH or HgCl2 (Hazama et al. 2002) AQP11 has been localized intracellularly in the proximal
tubule, and AQP11-null mice exhibited vacuolization and cyst formation in the proximal tubule,
demonstrating that AQP11 is essential for the function of this nephron segment (Morishita et al.
2005).
4
1.2
Aquaporin-4 (AQP4) in the Central Nervous System
Aquaporin-4 is one of the members of the aquaporin family, and was found to be
preferentially expressed in the central nervous system (CNS) (Hasegawa et al. 1994; Jung et al.
1994). AQP4 is the most abundant aquaporin in the astrocytes, especially in end-feet of the
processes which surround the brain capillaries (Frigeri et al. 1995; Nielsen et al. 1997; Nagelhus
et al. 1998; Venero et al. 2001; Tomás-Camardiel et al. 2004; Finnie et al. 2008). Astrocytes are
also the most abundant cell type in the CNS which protect and support normal neuronal activities
(Suzuki et al. 2006). Therefore, the location of AQP4 is specific to blood-tissue and tissue-CSF
border, possibly promoting fast water transport between those compartments (Bloch et al. 2006;
McAllister and Miller 2006; Tourdias et al. 2009). The expression level of AQP4 may influence
interstitial hydrocephalic edema and reducing the excess ventricular CSF (Mao et al. 2006;
McAllister and Miller 2006). In those literatures, the changes in AQP4 expression in specific
brain regions corresponded to the severity and duration of hydrocephalus.
1.2.1 Cytotoxic edema
In cytotoxic edema, the blood–brain barrier (BBB) remains intact. This edema is due to
the derangement in cellular metabolism resulting in inadequate functioning of the sodium and
potassium pump in the glial cell membrane. As a result there is cellular retention of sodium and
water. There are swollen astrocytes in gray and white matter. Cytotoxic edema is seen with
various intoxications. Verkman’s group disclosed that AQP4 null mice are protected from
several models of cytotoxic brain edema, including hyponatremia (Papadopoulos and Verkman,
2007).
5
Cytotoxic edema is caused by stroke and meningitis which is a condition of the lack of O2
supply to the brain. The lack of O2 reduces the amount of available ATP and prevents
transmembrane ion pumps (Na+) to buildup of Na ions in cell. Then it causes an influx of water
and increase the osmotic stress, which makes brain swelling and finally the malfunction of brain
tissues.
1.2.2 Vasogenic edema
Due to a breakdown of tight endothelial junctions which make up the BBB, vasogenic
edma occurs. It allows normally excluded intravascular proteins and fluid to perforate into
cerebral parenchymal extracellular space. Once plasma constituents cross the BBB, the edema
spreads very fast and widespread. The vasogenic edema is seen in response to trauma, tumors,
focal inflammation, late stages of cerebral ischemia and hypertensive encephalopathy.
The cytotoxic edema involves progressive cell swelling because of fast water uptakes,
whereas the vasogenic edema is the water release into the extracellular space due to the
malfunctions in the BBB. Mostly cytotoxic and vasogenic mechanisms coexist in the case of
water imbalance, typically one of both edema dominant in different brain desease; vasogenic
edema in tumors and cerebral abscesses, but cytotoxic edema in ischemic stroke and brain
trauma. (Zador et al. 2009)
1.3
A Morphology of Astrocyte Intracellular Space
Understanding cellular level structure is the essential springboard in implementing the
correct approach to model the intracellular behaviors. The cell structure includes subcellular
organelles such as a nucleus, endoplasmic reticular, golgi apparatus, etc. which are considered
6
compartments in a computational model and other dynamic molecules travel inward or outward
for the cell function. Each organelle or molecule is considered to play important roles in various
aspects of all possible cell functions. It is obvious that signaling and transporting molecules,
cytoskeleton-structures, and all existing organelles occupies a specific volume in the cell.
Approximately, up to 30% of the cellular volume is taken by the macromolecules (Ellis 2001)
and thus fundamentally the most prevalent transporting mechanism, diffusion, in a dilute
medium is expected to have a great hindrance from these crowding agents. The molecular
crowding is also considered to affect reaction rates of moving particle.
1.3.1 The Pathological Role of Astrocytes
Until quite recently, there has been a prevailing viewpoint on the astrocyte that those are
structural supportive glial cells in the brain and hence the dysfunction of astrocytes has not been
considered as a biological or pathological study. However current studies revealed that astrocytes
play a number of active roles in the brain, including metabolic support of neurons through
provision of glucose, the secretion or absorption of neural transmitters (astrocytes express
plasma membrane transporters for several neurotransmitters, including glutamate, ATP and
GABA), maintenance of the blood-brain barrier, and regulation of ion concentration in the
extracellular space (astrocytes express potassium channels at a high density) (Kolb and Whishaw,
2008; Araque et al. 1999). Moreover, astrocytes take part in the modulation of synaptic
transmission (through rapid changes in astrocyte morphology), vasomodulation, myelination
promotion (electrical activity in neurons causes them to release ATP, which causes astrocytes to
secrete cytokine leukemia inhibitory factor (LIF), a regulatory protein that promotes the
myelinating activity of oligodendrocytes), and nervous system repair (when injuries occur,
7
astrocytes become phagocytic to ingest the injured nerve cells, then fill up the space to form a
glial scar, repairing the area and replacing the CNS cells that cannot regenerate) (Piet et al. 2004;
Pascual et al. 2005). All these regulatory transmissions take places in their many processes.
In general, astrocytes are divided into two subtypes; fibrous and protoplasmic. Fibrous
astrocytes are found throughout all white matter and have relatively few organelles, and exhibit
morphology of long fiber-like processes. They often have "vascular feet" that physically connect
the cells to the outside of capillary wall when they are in close proximity to capillaries.
Protoplasmic astrocytes, found in grey matter in the brain exhibit morphology of several stem
branches that give rise to many finely branching processes in a uniform globoid
distribution. Both fibrous and protoplasmic astrocytes, when in proximity to the pia mater, send
out processes to form the pia-glial membrane (Sofroniew and Vinters 2010).
1.3.2 Water Homeostasis Functionality of Astrocyte
Numerous research since 1990 revealed that astrocyte processes are rich in the
aquaporin-4 (AQP4) water channel as well as K+, Ca2+, and N+ ions (Nedergaard et al. 2003;
Seifert et al. 2006). AQP4 proteins are densely polarized along astrocyte processes mainly
around the end-feet that contact blood vessels (Simard and Nedergaard 2004; Zador et al. 2009).
These proteins play a critical role in regulating water homeostasis in healthy CNS and play roles
in both vasogenic and cytotoxic edema as described in sections 1.2.1 and 1.2.2. Networks of
astrocytes linked together by gap junctions are thought to be able to rapidly dissipate small
molecules such as potassium and glutamate and prevent their potentially detrimental
accumulation. But this thesis limits on the AQP4 intracellular polarization into the plasma
membrane, the extracellular metabolism is beyond the scope and you can find more information
8
from somewhere else about AQP4 effect on extracellular metabolism (Kang and Othmer 2009;
Li et al. 2012; Nuriya et al. 2012; Wang et al. 2012).
1.3.3 The Organization of Microtubules
The rapid and directed transport of water as well as ions in the cellular level can be found
as microtubule (MT) based vesicle transport. Microtubules are filamentous intracellular
structures that are responsible for various kinds of movements in all eukaryotic cells. MTs are
constructed of long polymers with the diameter of ~25 nm and its volume fraction accounts for
about 3% of the cell volume (Luby-Phelps 2000).
Structurally, MTs consist of 13 protofilaments which are assemblies of tubulin in a headto-tail array. The tail ends (or minus ends) of MTs are usually anchored at MT organizing center
(MTOC). In astrocyte, most plus ends of MTs are extended to the processes to transport
materials to the end-feet as shown in Figure 1 and those positioning determines MT-dependent
vesicle transport.
At steady state, microtubules may appear to be completely stable, however there are
actions to shrink or grow constantly subject to regulation by cell as well as oscillation. For the
sake of brevity, only its structure was used as pathways of transporters and its functionality was
also limited to the unidirectional movement of kinesin motor proteins for the intracellular
transport as shown in section 1.5 and 2.3.
9
(a)
(b)
(a)
(c)
(b)
Figure 1. Organization of microtubules in various astrocytes; (a) immunofluorescence
image of microtubules of astrocytes using an anti-tunulin antibody (Tomás et al. 2005),
(b) Fluorescence speckle microscopy illuminating the intricate network f microtubules
(yellow) and actin filament (Purple) (image credit: Torsten Wittmann, UCSF).
1.4
AQP4 Synthesis Mechanisms
In general, protein synthesis is accomplished through a two-step process called DNA
transcription and translation. After a DNA is transcribed into a messenger RNA (mRNA)
molecule during transcription, the mRNA molecule is able to exit the nucleus through the
nuclear membrane pores and it enters right to the endoplasmic reticulum (ER) on which
ribosomes can be situated. Amino acids transferred by transfer RNA (tRNA) into the site of
protein synthesis are selected by opposing bases to the codons of the mRNA strand. This process
is called a translation. As a post-process of protein synthesis, the formed polypeptide chain
transported to golgi apparatus (GA) and experiences modifying, sorting, packaging into
macromolecule vesicles for the secretion into cytoplasm.
Water transfer-specific AQP4 protein synthesis is not different from the essential of
general protein synthesis procedure, which can be represented by simply two biological kinetic
10
reactions; transcription and translation. Transcription pathway is initiated by binding Nrf2 to
promoter regions containing the antioxidant-responsive element (ARE) motif in the aqp4 gene.
The Nrf2 bond gene turns to the active state for transcription of mRNA. As proteins are made up
of long chains of amino acids, AQP4 is encoded by the AQP4 gene by two different transcription
initiation sites (Jung et al. 1994). A longer isoform (34 kDa) exists in the brain and shorter
isoform (~31 kDa) can be found in other tissues such as lung. Umenishi and Verkman has
revealed how both isoforms are regulated when they are encoded from AQP4 cDNA (Umenishi
and Verkman 1998). In their work, the short isoform is encoded by the transcriptional event from
exon 1 to exon 4 as shown in Figure 2 and longer isoform is encoded by splicing event of exon 1
and transcriptional event from exon 0.
Figure 2. The organization of AQP4 gene and two isoforms of aquaporin-4 (Umenishi and
Verkman 1998).
1.4.1 Transcriptional Activation by Transcription Factor, Nrf2
Nrf2 is a transcription factor that regulates expression of many detoxification or
antioxidant enzymes. It is repressed by Keap1, Kelch-like-ECH-associated protein 1, which is a
cytoplasmic repressor that inhibits Nrf2 ability to translocate to the nucleus. Nrf2 and Keap1
interact with each other through the double glycine-rich domains of Keap1 and a hydrophilic
region in the Neh2 domain of Nrf2. Under quiescent condition, Nrf2 is constantly degraded via
11
the ubiquitin–proteasome pathway in a Keap1-dependent manner (McMahon et al. 2003;
Kobayashi et al. 2004). However, in response to oxidative stresses, Nrf2 degradation is ceased
and released from cytoplasmic Keap1 and translocated to the nucleus, and binds to antioxidant
response elements (AREs) in the promoters of its target genes to activate transcription. Although
the transcriptional regulators of AQP4 is not clearly found yet, Nrf2 is the one of the possible
regulators since the gene for NKCC1, like the AQP4 gene, is revealed to be affected by the
transcriptional activation of ARE (Cho et al. 2005).
1.4.2 AQP4 Transcription
The first part of AQP4 synthesis is transcription from DNA to mRNA. As well-known on
this procedure, a DNA sense strand is duplicated to mRNA strand from promoter sequence
where RNA polymerase attaches to when it reaches the terminator sequence. In this process,
AQP4 gene has two AQP4 transcripts; exon 0 promoter (M1) and exon 1 promoter (M23) (see
Figure 2). By different promoters’ locations, AQP4 mRNA is transcribed into two different
length isoforms. While Hirt et al. observed that M23 mRNA levels is higher than M1 mRNA in
mice brain, most quantitative immunoblotting and immunofluorescence studies tells M1 mRNA
is almost 8 folds higher than M23 mRNA and more increased in the ischemic condition
(Silberstein et al. 2004; Hirt et al. 2009).
1.4.3 AQP4 Translation
Two different AQP4 isoforms are encoded by their corresponding two different mRNA
isoforms. M1 isoform is encoded by M1 mRNA and M23 isoform is encoded by M23. Currently
new isoforms M23X (Zelenin et al. 2000) and M23A (Alikina et al. 2012) have been found, and
12
Holen and colleagues reported the existence of multiple isoforms in rat AQP4, including one that
formed a functional water channel, Mz ((Moe et al. 2008; Strand et al. 2009; Fenton et al. 2010).
It is known that Mz isoform expressed in rat but not human or mouse brain (Rossi et al. 2011).
But still it is considered that M1 and M23 are most dominant isoforms of AQP4 mRNA.
Ratio of AQP4 isoforms as well as translation efficacy in a cell are depending on the
expression levels of corresponding mRNAs. It is known that M1 AQP4 is encoded by one
mRNA, while M23 isoform can be encoded by several mRNAs. Even though there are several
factors affecting AQP4 gene expression, there is little information to explain which AQP4
mRNA isoforms are regulated.
1.4.4 Hydrocephalus Inducing AQP4 Up-regulation
Hydrocephalus is a pathophysiological abnormality associated with the increased amount
of cerebrospinal fluid (CSF) in a dilated ventricular system, i.e. cerebral edema. Since a CSF
consists mostly of water, and the AQP4 has been known to have a significant role in water
homeostasis in various systems as well as the central nervous system (Kimelberg 2004) since its
discovery, there have been numerous studies to reveal the close involvement of AQP4 in brain
edema. In detail, AQP4 is suggested that it plays a critical role in the early stage of brain water
accumulation caused by cytotoxic brain edema which is a type of the cerebral edema (Skjolding
et al. 2010). Solenov et al. reported that water permeability was seven-fold reduced in primary
astrocyte cultures from AQP4-deficient mice (Solenov et al. 2004). This result indicates AQP4
provides the predominant pathway for water movement in astrocytes.
However, unfortunately there exist only few clues that could help us understand the
pathophysiology underlying AQP4 up-regulation in hydrocephalus models. These studies tried to
13
indicate the association between aquaporin-4 and hydrocephalus in animal models like AQP4 knockout mice, kaolin- or LPC induced hydrocephalus in rats, and congenitally hydrocephalic
rats (H-Tx rats). Paul et al. evaluated increased expression of AQP4 in the congenital
hydrocephalus rat with advancing hydrocephalus (Paul et al. 2009). Shen et al. also studied the
expression of AQP4 in this animal model (Shen et al. 2006). AQP4 was highly expressed at the
ependymal lining, the end-feet processes of pericapillary astrocytes, and the subpial zone in HTx rats and indicated an adaptive mechanism. Shen et al. further discussed the role of AQP4 in
the pathophysiology of a subset H-Tx rat group with “arrested hydrocephalus”.
AQP4 was also up-regulated at the blood–CSF and blood–brain barrier interfaces in
hydrocephalic rat brains compared with control (Tourdias et al. 2009). Magnetic resonance
studies in these rats revealed a significantly bigger apparent diffusion coefficient (APC) and
larger CSF volumes, which were correlated with elevated expression of AQP4 under the
communicating inflammatory hydrocephalus.
Current primary treatment for hydrocephalus involves surgical drainage and diversion of
excess CSF. Verkman et al. showed that AQP4 deletion in mice accelerates the progression of
obstructive hydrocephalus which is caused by obstruction within the ventricular system, such as
a tumor, that prevents CSF proximal to the obstruction from draining into the subarachnoid space
where it is absorbed (Verkman et al. 2006).
1.4.5 Sulforaphane Enhancing Aquaporin-4 Expression
The vasogenic edema is a type of cerebral edema as a result of compromised bloodbrain barrier, thus induces fluid builds up in the intracellular space of the brain. Following a
traumatic brain injury (TBI), a decreased AQP4 is supposed to postpone the clearance of the
14
excess water out of the brain (Ke et al. 2001; Kiening et al. 2002). To confirm the AQP4
functionality under the brain edema disease resulting in the increased CSF volume in either
vasogenic and cytotoxic edema, Zhao et al. examined the AQP4 level in the injured brain (Zhao
et al. 2005). They found TBI decreased AQP4 level in the contusion core and post-injury
administration of sulforaphane (SFN), an isothiocyanate obtained from cruciferous vegetables,
attenuated AQP4 loss and further increased AQP4 protein levels in the penumbra region
compared with injured animals receiving vehicle. Up-regulated AQP4 in 3 days reduced the
water content to the level of 80% of it before the SFN injection.
1.5
Aquaporin-4 Transport Mechanisms
Broadly speaking, there are two main intracellular transport mechanisms; passive
diffusion in the cytoplasm of the cell and active motor-driven transport along microtubules.
Passive transport can be described in terms of the Brownian motion of free particles which
results in low speed transport. On the other hand, active transport is formulated in terms of
unidirectional motor protein motion along the cellular filaments for the long distance transport
with fast speed.
Diffusive motion of molecules still dominates the localization of proteins in the
cytoplasm although its motility is relatively lower than active transport. For normal diffusion, the
mean squared displacement of the molecules increases linearly with time, t, and the diffusion
coefficient, D, (1–20 µm2/s), and also depends on the dimensionality, d, (Elowitz et al. 1999;
Trinh et al. 2000; Van Kampen 2011). Elowitz et al. measured the diffusion of proteins in the
cytoplasm of Escherichia coli and found that proteins move randomly with a diffusion
coefficient of 3–8 μm2/s (Elowitz et al. 1999).
15
Shuttling protein molecules along the microtubule tracks is another major intracellular
pathway. Motor proteins (also called cargo proteins) have high affinity to the protein and lead to
convective fluxes that help rearrange the direct spatial localization of those molecules to the
plasma membrane (PM). Therefore this transport mechanism is fast and active. It is believed that
active transport is capable for the event of cellular damage, because the cell must answer
promptly by regulating relevant proteins to the targeted region but passive transport by diffusion
is not a prompt response.
Active transport along microtubules is managed by three different motor proteins; kinesin,
dynein, and myosin. Kinesin and dynein families move along microtubules whereas the myosin
family moves along actin filaments a shown in Figure 3. Kinesin is involved for the positive
directional movement which is toward the PM and dynein takes motion toward the negative
direction. Microtubules are polarized outwards to the peripheral of the cell membrane with their
‘plus’ ends from the center of the cell, with their ‘minus’ ends around the centrally located
microtubule organizing center (MTOC). The active transport of AQP4 proteins is a positive
directional movement from the center of the cell to the ‘plus’ ends of microtubules, so that
kinesin motor proteins and its convectional velocity must be considered.
16
Figure 3. Schematic representations of active bidirectional transport by motor proteins
along microtubules and active filaments in a cell. Myosin motor proteins walk along
actin filaments (red) at the cortex. Microtubule-based motors including the kinesin and
dynein motor proteins walk to the plus ends of microtubules (green), which are
oriented toward the actin cortex and to the minus end of the microtubule, which is
located at the microtubule-organizing center (MTOC, green) near the cell nucleus
(blue) respectively. Figure credit: (Ross et al. 2008)
17
2
METHODOLOGIES ON STOCHASTIC SIMULATION OF AQUAPORIN-4
INTRACELLULAR TRANSPORT
Most biological systems continuously exerting dynamic functionality take place at a
various range of spatial and temporal scales. Within the intracellular domain, its spatial scale
drops down to micrometers or nanometers, and same phenomena for the temporal scale. In this
re-confined domain, the number of molecules becomes low and events of molecular chemical
reaction take places with mille-seconds to minutes. Therefore, the principle continuum
approximation becomes invalid and the accumulated fluctuation from the low population may
bring about big off-track from the continuum model expectation. In other words, stochastic
effects in the micro-region become more important.
To fetch up the mathematical model with stochastic effects during aquaporin-4 (AQP4)
synthesis and polarization on the cell membrane needs two things; experimental data of AQP4
concentration on the specified region of synthesis and transportation or related prior knowledge
of the AQP4, aqp4 mRNA, and transcription factor such as kinetic reaction order, half-life,
concentration at steady-state in the intracellular systems and deterministic model without
stochastic effects. Most biochemical and molecular biology experimental data have a limited
accuracy due to errors, noise, incomplete information, and poor experimental design, which
contaminate the confidence of experimental data. Therefore, many hypotheses on the kinetic
modeling method to fit a proposed model to experimental data will be inevitably introduced
qualitatively to describe the biological systems and tested quantitatively to minimize errors
between a model and data. The stochastic model can be set up with two different ways by adding
18
stochastic terms into a validated deterministic differential system, which is called a general
stochastic differential system. Or random sampling using probability information of the events
from current state information and its kinetic model can realize the stochastic simulation.
In this thesis, we propose a model to investigate the synthesis of aquaporin-4 (AQP4)
protein and its transport mechanism to the targeting location, mainly end-feet of the astrocyte
processes. The synthesis pathway is limited with simple kinetic network controlled by feedback
transcription factor, Nrf2. The secreted AQP4 transport to the targeted area are assumed to have
two main transport phenomena; i) passive transport by free diffusion regardless of concentration
bias, and ii) active transport along microtubules. The compartments consisting the intracellular
environment in the model are endoplasmic reticular (ER), and nucleus, and microtubules (MT).
In our model, we also made the following assumptions. (i) Ribosomes, RNases, and other
factors involved in protein degradation and synthesis have long half-lives relative to the mRNA
or AQP4, and their effects on AQP4 regulation are omitted. (ii) AQP4 translation is constitutive,
and its degradation rate does not vary significantly over the cell. (iii) Since the mRNA transcripts
serve as a point source for the proteins for which they code, there is no spatial transport of AQP4
until a AQP4 vesicle made up with numbers of AQP4 molecules leaves ER membrane. (v) The
process of translation is such that the probability of proteins being made in the time interval (t,
t+dt) is independent of the number of proteins made before time t. (vi) The process of
degradation is such that the probability that a protein is degraded in the time interval (t, t+dt) is
independent of the number of proteins degraded before time t.
19
2.1
Intracellular Structure and Conditions
In the present study our aim is to obtain a mathematical model of intercellular
polarization of AQP4 in an astrocyte cell and in this dissertation, and therefore we are focusing
on the compartments and molecules which have relatively important functionality on AQP4
synthesis and trafficking activities. The compartments include a nucleus, endoplasmic reticula
(ER), Golgi Apparatus (GA), cell plasma membrane (PM), cytoskeleton, and most spatialdominant cytoplasm as shown in Figure 4. These structural units are taken into the modeling as
compartments and different stations for the generation of AQP4 from nucleus and its trafficking
in cytoplasm to plasma membrane (PM). The average diameter of a nucleus is approximately
6μm, which occupies about 10% of the total cell volume (Bruce Alberts 2002).
Figure 4. Sketch of intracellular structure of astrocyte.
The ER represents the entry point for newly-made proteins destined for the secretory and
endocytotic pathways. The ER is a membranous labyrinth located throughout the cytoplasm of
the cell as a single continuous network of flattened sacs and tubules. The total volume fraction of
20
ER takes up to 10% of cell volume as large as the nucleus (Voeltz et al. 2002). The role of ER is
the modification, folding, and assembly of proteins before the vesicle formation. Until it is
incorporated into transport vesicles budding off from the ER, it will stay in the ER to be folding
and oliogomer assembly. The coated vesicles in ER go to GA by membranous transport. In GA it
is known that proteins become other bigger complex protein (Hong 1998). The ER-GA route is
considered transitional transport between ER and GA membranes. But, since mechanistic
involvement in ER-GA transport remains debatable, we could not model exact pathway from ER
to GA. Instead, in our proposed model, we lumped ER and GA unit into a singular spherical
compartment as a transitional element and placed a centriole on the boundary of this
compartment. Centrioles are known as the origin of microtubule structures in the cytoplasm
(Beisson and Wright 2003; Feldman et al. 2007).
Other important compartments are cytoskeletons throughout the cytoplasm. Main
functionality of cytoskeletons is known as structural supports or anchorage for intracellular
organelles. But current researches found cytoskeletons have a role of active and directed
transportation of molecules; for examples signaling molecules along the actin microfilament
from PM to the link with microtubules (Kerr et al. 2003) and secretory molecules along
microtubules to PM (Kholodenko 2003). The cytoskeleton consists of microtubules (diameter
~25 nm; length ~100 μm), intermediate filament (diameter ~10 nm), and actin filament (diameter
~7 nm; length 30-100 μm). The volume fraction of all cytoskeleton takes about 3% of entire cell
volume (Luby-Phelps 1999). Since the cell molecular simulator in this study only counts on the
intracellular transport of positive direction to PM by kinesin motor proteins, we implemented
microtubules (MTs) in our 3d model. The velocity of motility generated by kinesin-1 is reported
21
between 0.4 and 0.9 µm/s (Vale et al. 1985; von Massow et al. 1989; Steinberg and Schliwa 1996;
Böhm et al. 1997). The detailed information of MT modeling is described in section 2.1.1. All
compartments considered in the cell molecular simulator is listed in Table I.
Table I. Fixed compartments design specification
Compartments
dimension
Astrocyte cell diameter
~15 μm
Length of astrocyte process
~15 μm
Astrocyte cell nucleus diameter, dnucleus
~6 μm
Endoplasmic reticulum outer diameter
~1.5 μm
No. of cytoskeleton microtubules, Nm
3300
Length of microtubules, Lm
http://users.rcn.com/jkimball.ma.ultranet/BiologyPages/C/Cytoskeleton.html
Diameter of microtubules, dm
http://users.rcn.com/jkimball.ma.ultranet/BiologyPages/C/Cytoskeleton.html
Volume fraction of microtubules, Vm
(calculated)
~25 μm
~25 nm
1.15%
2.1.1 Aquaporin-4 Water Channels in Cell Membrane
Aquaporin-4 water channel in the cell membrane is referred to orthogonally arranged
particles (OAP) in which two different isoforms of aquaporin-4 proteins, M1 and M23, compose
of one orthogonal supramolecular tetramer. In principle, AQP4-M23 is the OAP-forming
isoform since the projection from M23 monomer is a potential site of M23-M23 intertetrameric
interaction to build large OAP, whereas AQP4-M1 alone is unable to form OAPs due to the lack
of bridges for M1-M1 interaction but can coassemble with M23 in OAPs as heterotetramers.
Therefore the modeling of aquaporin-4 polarization forming 6 different possible OAP water
22
channel from M1 and M23 aquporin 4 monomers must consider at least these two different
monomers as moving particles and the expression ratio of water protein channels to both AQP4
monomers should satisfy 1:4. The probability of M23 synthesis was set 3 fold more abundant
and there is no differences of M23-M23 and M23-M1 interaction probability (Jung et al. 1994;
Lu et al. 1996). The experimental quantitative analysis was driven by quantum dot single particle
tracking, single molecule photobleaching, and blue-native gel electrophoresis (Crane et al. 2009;
Tajima et al. 2010).
The ratio of M1/M23 composition was determined using random numbers (for each
monomer). Newly generated AQP4 tetramers were then appended to an existing OAP until no
more tetramers could be added. The rigorous mathematical model of OAP assembly was
explained by (Jin et al. 2011) and the schematic illustration is shown in Figure 5.
Figure 5. Schematic of OAP assembly from AQP4 M1/M23 heterotetramers. (A) (top) M23 and
M1 AQP4 monomers. Projection from M23 monomer denotes potential site of M23-M23
intertetrameric interaction. (bottom) Six possible independent AQP4 M1/M23 heterotetramers.
Four-way rotational variants of each are possible as well. (B) Examples of OAPs at low and high
M23:M1 ratios. (Jin et al. 2011)
23
Individual AQP4 tetramers contain zero, one, two, three, or four M1 monomers, in six
possible independent configurations (with rotation possible) (Figure 5.A). If tetramer-tetramer
association requires at least one intertetrameric M23-M23 interaction, then many OAP
configurations are possible of different size, shape, and composition. Conceptually, more M1 at
fixed M23 is predicted to reduce OAP size because M1 blocks tetramer associations by coating
the outer surface of OAPs (Figure 5.B). The modeling here involved the serial generation of
AQP4 tetramers whose M1/M23 composition was determined using random numbers (for each
monomer) and the specified M1:M23 ratio.
<Assumption>
From the literature, we have known that each AQP4 will have three different secondary
structures by attractive and repulsive force between amino acids and the expression ratio of
heterotetramers of OAPs is species-dependently different. But in this modeling, we assume the
secondary structures of all AQP4s are varying in terms of their size when they form a vesicle
because different composition itself does not have a significantly different functionality in the
formation of OAP and water transport performance. Due to the species-dependency of OAP, the
generic model for its localization is not possible and this generalization for human model is
beyond of this thesis.
2.1.2 Modeling of Microtubule Structure
Microtubules (MT) exhibit dynamic instability with growth, shortening, and movement.
Furthermore, assembled tubulin structures are known to have various conformations. Therefore
to present MT structures more accurately, it needs the intensive study for the assembly and
24
disassembly mechanisms and it is beyond of the scope in this thesis. To represent MTs’ structure
and incorporate with 3d astrocyte model, we modeled the static structure of microtubules. For
this task, we developed smooth & shortest path finding algorithm, which aims a finding the
smooth and shortest path in the presence of obstacles.
This algorithm begins with a random directional growth of MT segment from MTOC to
the neighboring mesh. From 2 nd growth, it uses the directional information from target vector,
A(t), which is a unit vector connecting a current point and selected end point. To keep the
smooth curvature, it also takes the information of tangent vector, T(t). The next point by each
growth step will be determined by the summation vector, D(t), of both vectors with small
variation by random vector.
Pi ( x, y, z )  Pi 1 ( x, y, z )
Pi ( x, y, z )  Pi 1 ( x, y, z )
(1)
Ptarget ( x, y, z )  Pi 1 ( x, y, z )
(2)
Ti  a 
Ai  b 
Ptarget ( x, y, z )  Pi 1 ( x, y, z )
Once we have a vector, D(t), its position will have small variation to avoid an obstacle or a
breakaway out of cell membrane.
During the growth of MTs, obstacles are given as region with the range of dimension.
When it reaches the region of obstacle, the algorithm finds alternate directional segment
increasing random vector effect and relatively decreasing the effect of the summation vector of
both tangent and target vectors. For example, Figure 6 shows 2d smooth & shortest path finding
algorithm test. In that figure, all growths started from the center and made a detour from four
rectangle shape obstacle regions and arrived to the four selected target point.
25
Figure 6. An example of Smooth and Shortest Path Finding algorithm in 2d. All paths avoid the
obstacles (rectangular shape) on the straight way to the target points (marked as ‘x’).
This algorithm also prohibits the pathways to breakaway using normal vector and
distance information to the closest boundary of cell. When it approaches to the boundary from
the distance information, the increase of normal vector effect will offset the momentum to the
boundary and gradually turns over the direction of growth to the target point. The procedure of
this algorithm is described as follows.
Smooth & Shortest path finding algorithm

Step 1. Select a starting mesh of MT growth (assumed as centriole) and target mesh
positioned at the end-feet.

Step 2. Generate a first growth to random direction with the length to the connecting
mesh.
26

Step 3. Grow a segment to the neighboring mesh adjusting its direction and length
with the information of tangent vector and target vector (direction to the final
destination)

Step 4. Repeat step 3 until it does not meet the boundary of process.

Step 5. If the growth does not enter to the process boundary, the growth will stop and
new MT growth will start from step 1. Otherwise, the segment of MT goes into the
process which contains a target mesh selected in step 1.

Step 6. Grow a segment to the neighboring mesh adjusting its direction and length
with the information of tangent vector, target vector, and distance to the normal
surface.

Step 7. Repeat step 6 until it comes to the end-foot boundary.

Step 8. When it comes to the mesh which is one of end-feet members, the smooth
path finding algorithm will stop and finalize its growing to the target mesh selected in
step 1.

Step 9. The end criterion for the total number of MTs will terminate the algorithm.
Note 1. Each step of growth can record its position as mesh number or absolute position.
Note 2. Each step of growth will have slight directional and longitudinal variation to
avoid when it approaches to an obstacle or boundary of cell. The combinatory consideration of
tangent vector for smoothness and distance to the closest surface of current point for shorter
pathway will gradually modify the direction of growth.
27
2.2
Stochastic Model and Simulation of AQP4 Synthesis
2.2.1 Kinetic Reaction Model of AQP4 Synthesis
The synthesis of AQP4 water channels roughly consists of the transcription and
translation mechanisms. While entire mechanism of AQP4 synthesis is such a complex reaction
system, we used a simplified lumping model. The lumped model does not cope with the details
of innumerable components but provides the reasonable information on factors we focus on the
model to regulate AQP4 proteins in cellular level. Since we are interested in the AQP4 regulation
corresponding to the change of mRNA molecule by the external stimulus, other factors are
disregarded in this thesis.
The entire kinetic reaction mechanism of AQP4 synthesis in an astrocyte can be initiated
with a suloforaphane (SFN) administration. SFN stimulates transcriptional activation through the
binding of Nrf2 to the promoter regions containing the antioxidant-responsive element (ARE)
motif of the aqp4 gene. The transcription of mRNA transcripts occurs with a higher frequency
(when Nrf2 is bound) than the basal transcription rate of the aqp4 gene. In our mathematical
formulations, the Nrf2 bond gene turns to the active state for transcription of mRNA.
The transcription rate of aqp4 mRNA is determined with ratio of transcription from gene
and level of activated gene. The SFN treatment increases Nrf2 release into nucleus and binds on
aqp4 gene promoter to shift transcription rate fast. Until 8 hours after SFN treatment, upregulated transcription continues and its rate constant returns to the initial steady state level with
Nrf2 exportation to cytoplasm after 8 hours. In this model, translation kinetics is 1 st order
reaction of mRNA concentration and all degradation kinetics are set to follow 1 st order reaction.
The overall lumped model is depicted in Figure 7.
28
SFN: sulforaphane
DNA*: activated DNA
Nrf2: transcription factor
AQP4: aquaporin 4 protein
Keap1: cytoplasmic repressor of Nrf2
Keap1
Nrf2
SFN
Nrf2 Localization
DNA Activation
DNA
DNA*
Nrf2
DNA Inactivation
ø
mRNA transcription
mRNA
AQP4 translation
Degradation
ø
AQP4
Figure 7. Schematic mechanism of AQP4 protein synthesis. In this study, Nrf2
activation by SFN treatment up-regulates AQP4 expression on the astrocyte cell.
The procedure of protein synthesis and its transport pathways are not reproducible due to
the stochastic variations (Arkin et al. 1998). Therefore, a continuous simulation will not be
appropriate for the current system and the proposed cellular simulator used stochastic model and
simulation to capture discrete events of reaction and transport mechanisms. In general, the
probabilistic representation with chemical master equation (CME) has been applied to many
cases such that gene transcription, intracellular signaling pathway, molecular diffusion and
reaction, and etc. (Baras and Mansour 1996; Paulsson et al. 2000; Khanin and Higham 2008;
Isaacson and Isaacson 2009). But the CME method is often very hard or impossible to solve
analytically and numerically due to the rapid increase of computational duty when the
dimensionality of system goes more complex since the CME account for all possible tracks of
events (see section 2.3.2). Therefore we used an approximation of CME method using Gillespie
29
algorithm and time discretization for stochastic simulation. The entire simulator process is
summarized in the Figure 8.
Figure 8. Cell molecular simulator flowchart.
30
In reaction mechanisms, we incorporated propensity functions to calculate integer
number of reaction events, in which propensity functions, a, is updated each simulation time
based on the current state, j, population of species as shown in eqs. (3) to (5) for the synthesis of
AQP4 and eqs. (6) to (8) for the inactivation or degradation of the species. The discrete numbers
of events for each reaction are then sampled by Poisson distribution. To avoid over-estimating
the total number of events by Poisson distribution, the discrete simulation step size needs to be at
least 10 times lower than the minimum rate propensity among all possible reaction events.
a( j , DNA* )  kb [ DNA] j [TF ] j
(3)
a( j , mRNA)  ktc [ DNA* ] j
(4)
a( j , AQP4)  ktl [mRNA] j
(5)
a( j , DNA)  kub [ DNA* ] j
(6)
a( j , mRNA0 )  kdr [mRNA] j
(7)
a( j , mRNA0 )  kdr [mRNA] j
(8)
The stochastic model and simulation on the synthesis of AQP4 is also verified with a
deterministic model by averaging repeated stochastic simulations to offset the stochastic
fluctuations. The mathematical model with ODEs for the biological kinetic ordinary differential
equations (ODEs) is described in eqs. (9) to (12).
d [ DNA]
 kb [ DNA][TF ]  kub [ DNA* ]
dt
d [ DNA* ]
 kub [ DNA* ]  kb [ DNA][TF ]
dt
(9)
(10)
31
d [mRNA]
 ktc [ DNA* ]  kdr [mRNA]
dt
d [ AQP 4]
 ktl [mRNA]  kda [ AQP4]
dt
(11)
(12)
2.2.2 Transcription Kinetics
Transcription pathway is initiated by binding Nrf2 to promoter regions containing the
antioxidant-responsive element (ARE) motif in the aqp4 gene. The Nrf2 bond gene turns to the
active state for transcription of mRNA. The transcription rate of aqp4 mRNA is determined with
ratio of transcription from gene and level of activated gene (eqs. (14) and (15)). The SFN
treatment increases Nrf2 release into nucleus and binds on aqp4 gene promoter to shift
transcription rate fast. Until 8 hours after SFN treatment, up-regulated transcription continues
and its rate constant returns to the initial steady state level with Nrf2 exportation to cytoplasm
after 8 hours.
d [ DNA]
 kb [ DNA][TF ]  kub [ DNA* ]
dt
(13)
d [ DNA* ]
 kub [ DNA* ]  kb [ DNA][TF ]
dt
(14)
d [mRNA]
 ktc [ DNA* ]  kdr [mRNA]
dt
(15)
d [ AQP4]
 ktl [mRNA]  kda [ AQP4]
dt
(16)
Schwanhausser et al. quantified the transcription rate of two mRNA molecules per hour
as a median rate (Schwanhausser et al. 2011), and other microscopic study on the
32
cytomegalovirus (CMV) promoter reported transcription termination rates of 5.8 to 8.7 mRNAs
per hour (Darzacq et al. 2007).
2.2.3 Translation Kinetics
mRNA transcribed from the gene is ready to bind of ribosome to the ribosomal binding
site and then ribosome facilitate the assembly of the polypeptide chain for a single protein or
multiple proteins at a same mRNA macromolecule simultaneously. The entire process from
initiation with ribosome to termination with released complete peptide is called protein
translation. The advanced mathematical frameworks on this process have been proposed by
many researchers (MacDonald et al. 1968; Lee et al. 2003; Mehra and Hatzimanikatis 2006), but
the characteristics of codons, tRNAs, elongation factors, and etc. as well as the intermediate
complex steps of elongation cycle kinetics respect to the ribosomal state are beyond of this
thesis’s scope. Therefore, as similar to transcription kinetics in previous chapter, our
mathematical model will be derived from simple principles of first-order AQP4 translation and
its first-order decay.
2.2.4 Degradation Kinetics
Degradation of mRNA and AQP4 is an important mechanism in determining the levels
and regulation of gene expression. As a feedback control pathway of forward procedure, it
understandably is responsible for the actual regulation of corresponding molecules in the
network. For example, the kinetic behavior of individual RNA polymerase II (RNAPII)
transcribing a gene provides a precise quantification of the contribution of mRNA synthesis to
the cellular pool of transcripts (Trcek et al. 2011). Generally, modified cells inhibiting any
33
possible transcription processes is used to quantify a mRNA decay and obtain kinetic
information, but the detailed mechanism of degradation is beyond scope of this dissertation.
Instead, we adopt first-order kinetics to explain mRNA and AQP4 degradation with 4.8 hours
and 24 hours half-life respectively. Simply degradation rate constant can be calculated with
following relations with half-life in the case of first-order reaction.
kdegradation 
2.3
0.693
t1 2
(17)
Stochastic Simulation of the Intracellular Transports
Most biological systems continuously exerting dynamic functionality take place at a
various range of spatial and temporal scales. Within the intracellular domain, its spatial scale
drops down to micrometers or nanometers, and same phenomena for the temporal scale. In this
re-confined domain, the number of molecules becomes low and events of molecular chemical
reaction take places with milliseconds to minutes. Therefore, the principle continuum
approximation becomes invalid and the accumulated fluctuation from the low population may
bring about big off-track from the continuum model expectation. In other words, stochastic
effects become more important. Looking into last two decades, in the field of molecular biology,
the stochasticity of the system and processes occurred in the system has been increasingly
considered and it led to a better understanding of biological systems at the intracellular scale.
Therefore, the intracellular kinetics between molecules such as DNA, RNA, and protein should
consider additionally beside an important law of physics; randomness and discreteness.
Intracellular processes with lower number of particles such as extreme case of single DNA
template in the gene transcription process could not be interpreted with ratio of probability so
34
that discrete events on that particle with randomness caused by probability may be explained in a
given time scale. This discreteness of molecular events also influences the fluctuation of the
number of molecules (or concentration) and accordingly causes intrinsic noise of the biochemical
interactions.
In principle, the dynamics of intercellular transport of AQP4 protein can be analyzed with
ordinary differential equations (ODEs) and can be solved with continuous simulation in the
temporal integration. In order to analyze transport mechanism through the cell, a particle-based
simulation has been developed, which preserves the stochastic nature of reactions as well as the
undirected motion of diffusion alongside microtubule tracks. The spatial location of each
molecule is tracked in the simulation through a true cell model. As the building of true cell and
transport mechanistic model, a couple of scenarios of aquaporin-4 regulation is simulated and
validated with the experimental expression of aquaporin-4 in the astrocyte end-feet. This
regulation is future researched with the hydrocephalic brain disease.
In this section, I focused on the spatiotemporal distribution of the aquaporin-4 molecules
and its stochastic probability of the concentration on the plasma membrane as the final
destination of transport. The intracellular transport pathway will be identified and applied in the
stochastic simulation to determine the probability of event occurred in intracellular
compartments such as ER, GA, and cytoplasm.
2.3.1 Mechanisms of AQP4 Transport
The molecules secreted into cytoplasm from nucleus will have two different transport
mechanisms by passive pure diffusion and active motor-protein driven convective transport. For
the stochastic simulation of pure diffusion, we used fractional Gillespie multi-particle algorithm
35
(fGMP) which is an approximation of Gillespie method with discrete time interval. This method
is also based on the Gillespie multi-particle (GMP) method proposed by Rodriguez et al.
(Rodríguez et al. 2006), which is a particle-based spatial stochastic method to simulate
biochemical networks. In the spatial pure diffusion problem, fGMP method computes four
different directional diffusion group by normal distribution, N(μ,σ), with the distances to
neighboring sub-volumes. Assuming the homogeneous diffusivity, 0.125 μm2/s, all diffusion
probabilities are inverse proportional only to the distances.
The active transport by motor protein along microtubules has three mechanisms;
association to microtubules, dissociation from microtubules, and convective transport. The
motor-protein is highly attractive of the diffusing molecules and the dissociation of molecules is
relatively low. The normalized convective transport can be calculated from eq. (18) (Vale, 2003).
(18)
The length of the step depends on the velocity of the motor protein vmotor (~0.5 μm/s) and the
active transport is deterministically simulated with stochastic association and dissociation
reaction events (Vale et al. 1985; von Massow et al. 1989; Steinberg and Schliwa 1996).
For a rigorous simulation, entire transport network can be broken into single transport
between two compartments which can be named donor and receiver of vesicles. Each
compartment is assumed to have evenly spatial-distribution of tracking molecules so that updates
of stochastic event do not consider the internal movement of molecules. Simply, the transport
network will be modeled as shown in Figure 9.
36
End-feet
Dissociation
Reaction
fusion
rate
Receiver
vesicle
Dissociation
Reaction
active transport
by motor protein
Inactive transport by diffusion
Association
reaction
Budding rate
Donor
Figure 9. Simplified transporting cascade of protein molecule between compartments.
The procedure of vesicle transport from donor to receiver is described below.
Budding Process
Step 1. Donor compartment buds transporting molecules in the vesicles
Transporting Process
Step 2. Budding vesicle diffuses into cytoplasm.
Step 3. Budding vesicle transports by inactive diffusion.
Step 3-1. Vesicle associates to the motor protein.
Step 3-2. Vesicle associated to motor protein transports by active motor protein transport.
Step 3-3. Some vesicle associated to motor protein dissociates from motor protein during
the cytoskeleton track.
Step 3-4. Some vesicle associated to motor protein dissociates from motor protein at the
end of the cytoskeleton.
Step 4. Diffusive or transporting via motor protein arrives at receiving compartment.
Fusion Process
Step 5. Protein molecules fused from vesicle to receiving compartment.
37
On the first look, vesicle budding looks like a simple reaction in which new vesicles are
created with a specific rate. However, this rate is determined by the respective cargo for the new
vesicle (on demand) and furthermore by the availability of the respective machinery molecules.
In the simulation, we assumed that vesicle size varies from 1 to 50 molecules and budding
process was initiated randomly since there is no information of formation rate of a vesicle - cargo
dimer.
Fusion process can be modeled as a simple reaction between two agents. In spatial
modeling, the distance of reacting agents should be placed in a reasonably small grid and the
reaction rate determines the probability of reaction. In reality, the fusion process terminates the
intracellular transport of AQP4 molecules in this study. So we set very low reaction activity
when the vesicle is not close to end-feet and relatively high reaction activity when it reaches the
end-feet of the processes.
2.3.2 The Chemical Master Equation
The most typical method to describe the random events is probabilistic representation
with chemical master equation (CME). CME has been applied to many cases such that gene
transcription, intracellular signaling pathway, molecular diffusion and reaction, and etc. (Baras
and Mansour 1996; Paulsson et al. 2000; Khanin and Higham 2008; Isaacson and Isaacson
2009). To make the description of CME clearer, we included an example of reaction system.
Consider the following system of reactions.
k1
A 
B
k2
B 
C
(19)
38
with rate constants k1 and k2 given as 0.01 and 0.005 respectively. The macroscopic description
of how the amount of A, B, and C changes with time would be expressed with the evolution
equations as
dA
  k1 A(t )
dt
dB
 k1 A(t )  k2 B (t )
dt
dC
 k2 B (t )
dt
(20)
Note that an approximation is being made since we vary concentrations continuously with time.
In this system, with time interval ∆t, the probability function at time t+∆t can be obtained as
eq. (21).
m
dP( X, t  dt )
 P( X, t ) P(no change over dt )   P( X  v r , t )P( state change over dt )
dt
r 1
(21)
where, vr is a stoichiometric vector and X is a concentration vector. The probabilities of ‘no
change over ∆t and state change over ∆t are expressed as eqs. (22) and (23).
P(no change over dt )  P(no reaction over dt )
m
 1   ar ( X)dt
(22)
r 1
P( state change over dt )  P(reaction r occurs )
 ar ( X  v r )dt
(23)
Eq. (21) can be re-written as eq. (24).
m

 m
P( X, t ) 1   ar ( X)dt    P( X  v r , t )  ar ( X  v r )dt   P( X, t )
P( X, t )
 r 1
 r 1

t
dt
m
m
  P( X, t ) ar ( X)   P( X  v r , t )ar ( X  v r )
m
r 1
r 1
   P( X  v r , t )ar ( X  v r )  P( X, t ) ar ( X) 
r 1
(24)
39
2.3.3 Monte Carlo simulation and the Gillespie algorithm
The CME method in the previous section 2.3.2 is often very hard or impossible to solve
analytically and numerically. Solving CME encounters a rapid increase of computational duty
when the dimensionality of system goes more complex since the CME account for all possible
tracks of events. Therefore many approximations of CME model have been proposed to avoid
unnecessary computation of probabilities by selecting one series of events in a single scenario
from all possible events. Gillespie’s algorithm is one of most famous method to simulate a
realization (trajectory) of the system with correct statistics (Gillespie 1977). In fact, Gillespie
algorithm is a good approximation of CME using uniformly distributed randomness of 0 or 1
event for the all possible reaction events. Gillespie algorithm substantially record a single event
based on the current state of system, and continuously tracks a series of all single events.
Therefore, Gillespie algorithm still needs high duty of computation to track the scale of millions
of intracellular events.
To trade-off between computational duty and accuracy, other approximations are
proposed with leaping of intracellular time interval such as tau-leaping (Gillespie 2001), rleaping (Auger et al. 2006), k-leaping methods (Cai and Xu 2007), and next reaction method
(Gibson and Bruck 2000), and Gillespie multi-particle method. Among these approximated
methods, we selected Gillespie multi-particle method since it is based on the Diffusion-Reaction
Master Equation (RDME), and also capable to simulate diffusing particles either one by one or in
groups while conserving the normal diffusion properties. In general, diffusion events in a cell
take place relatively so often than other reaction events since diffusion is not limited from the
distance of reacting molecules and the affinity between each other. Therefore grouping diffusive
40
particles can drastically reduce the computational cost as the dimensionality of molecular
population goes higher.
2.3.4 Gillespie Multi-Particle (GMP) Method with Modified Diffusion Fraction for the
Structured Coordinate System
To approximate the stochastic diffusion events in this work, we used Gillespie multiparticle (GMP) method proposed by Rodriguez et al. (Rodríguez et al. 2006). This is particlebased spatial stochastic method to simulate biochemical networks. It does not compute full
pathways which demand computational expensiveness and is possible to deal with spatial effects
that ODE model could not cope with. This method works for the heterogeneous distributed
system by discretizing space into smaller homogeneous sub-volumes. Then homogeneous subvolumes then take stochastic events of reaction or diffusion updating current sub-volume states.
In this study, we are proposing a mathematical model to predict AQP4 protein polarization on
the astrocyte end-feet, and it needs to direct comparison of state corresponding to the experiment
analysis, so the generic GMP method is modified with discrete time interval.
In the spatial diffusion problem, GMP method computes four different directional
diffusion group by normal distribution, N(μ,σ) as shown in for the case of Cartesian structured
coordinate system. Assuming the homogeneous diffusion in the structured mesh, all diffusion
probabilities becomes identical as *Pn=Ps=Pe=Pw=0.25 (*Pn : probability of diffusion to north
direction).
41
Figure 10. Structured meshes in the Cartesian coordinate system for different
directional diffusion
Therefore average times to travel to neighboring sites are not different but the randomness of
event and discreteness of molecules make the fluctuation of number of molecules in the
neighboring sites after specified time interval (∆τ) and its conceptual expression is depicted in
Figure 11. This method computes approximated number of diffusing molecules to the
neighboring sub-volumes at a single simulation time with full clearance of origin sub-volume.
This empty sub-volume can be avoided by incoming molecules diffused from neighboring subvolume(s), but initial diffusion event from source point or the distributed system where there is
small number of molecules may not avoid this undesired empty site problem.
Figure 11. Diffusion event simulation by Gillespie multi-particle method (GMP). The
fluctuation of diffusing molecules is caused by randomness of event and discreteness
of number of molecule. In the original GMP method does not account for the empty
site situation after single diffusion event simulation when a group of molecules is
initially located only in the specific sub-volumes or where a site has all empty
neighboring sites.
42
Therefore, in this study we proposed modified fractional GMP method which computes
some fraction of diffusion with shorter event interval (∆τ). This method is very effective in the
structured system as shown in Figure 12.
Figure 12. Diffusion event simulation by modified fractional Gillespie multi-particle
method. Fractional diffusion prevents the undesired empty site problem happened at
initial steps of diffusion by source point or when a site are surrounded with all empty
neighboring sites. Modifying diffusion ratio from 1 to 0.8 (typically we use 80% of
clearance of diffusion event site) shortens simulation interval from ∆τ to 0.8*∆τ.
2.3.5 Gillespie Multi-Particle (GMP) Method Event Time Computation for the
Unstructured Coordinate System
We have known that discretized time interval in the structured coordinate system because
the traveling time to the neighboring sites for all sub-volumes are identical. However, to describe
the biological system such as astrocyte in this study, we may necessarily use unstructured
meshes due to the computational purpose or geometric factors. Therefore each sub-volume has
different traveling time (eq. (13)) due to the different geometric distance to neighboring subvolumes. The event times for all sub-volumes are initially computed and simulation follows the
increasing order of event time. Then next event time for a sub-volume where an event fired is
updated with the interval of initial event time.
43
xy
2dDs
d : system dimensionality, Ds: diffusion coeff. of species S
s 
(25)
When we consider stochastic diffusion-reaction simulation, the list of event times of all
sub-volumes is initially calculated based on the diffusion events from average travelling time,
and then the combinatorial simulation between diffusion and reaction is performed with
upcoming next event. This process is explained with conceptual flowsheet in Figure 13.
Setup meshes and
initial particle distribution
Initialize all species diffusion time steps
Calculate first diffusion time, τD=min{τi}
Calculate first reaction time τR
τD < τR ?
Record tsim=τD
Record tsim=τR
Diffuse species i
for all sites at τD
μ Reaction fires
at time τR
Update τi=τi+τi
Update tsim=tsim+τR
Next diffusion time
τD=min{τi}
Next reaction time
τR
no
τR > tend ?
τD > tend ?
yes
no
End
Figure 13. Stochastic Reaction-Diffusion Flowsheet. The fractional diffusion is not considered in
this conceptual flowsheet.
44
3
RESULT OF THE AQP4 INTRACELLULAR TRANSPORT BY STOCHASTIC
SPATIOTEMPORAL REACTION-DIFFUSION SIMULATION
3.1
Astrocyte Cell Reconstruction
The rationale of mathematical modeling of intracellular transport is to integrate
experimental observations of cell morphology (confocal microscopy), transcription factor
activation (immunofluorescence), and induced target protein expression (western blot) using a
cellular simulator based on stochastic simulations. Therefore the cell model reconstruction from
cell morphology taken by clinical image processing is necessarily a first task of the cell
molecular simulator. Using mimics mesh generator (Materialise), the geometry of the cellular
simulator was reconstructed from a stack of confocal images of a single astrocyte. Figure 14 (a)
shows that microscopic image of a single astrocyte and Figure 14 (b) is the result of
reconstruction to 3D model with 11,390 volume meshes.
Figure 14. Cell reconstruction. The series of confocal microscopic images (a) are reconstructed to
the computational model (b). The computational model has 11,390 tetrahedron volume meshes.
45
To reduce the noise in the microscopic images, we applied median filter. The
unnecessary regions which are not the parts of cell in the filtered image are then reduced and the
indistinct parts due to the resolution of microscopic images are masked to be surely added when
the 3D surface mesh is created. However, after creation of 3D volume meshes without any postprocessing, the models often do not meet the required quality. To improve the quality, we postprocessed the smoothing and wrapping with proper index of number of smoothing iteration and
wrapping distance in the mesh generation software. The improvement applied in the model is
visualized in Figure 15.
Figure 15. Cell meshing post-processing of wrapping and smoothing. Selecting proper
parameters in Mimics mesh generator, we can obtain the high quality 3d model of astrocyte
cell.
3.2
Microtubule Compartment Implementation
The confocal microscopic images used in cell modeling do not show all necessary
compartments such as nucleus, ER, and cytoskeletons. Those compartments are important to
place budding molecules from ER membrane and set pathways of active transport. To simulate
46
natural transport mechanisms, we implemented invisible nucleus, ER and microtubules. For the
nucleus and ER, we used CAD mask option in Mimics software to add those regions into 3D
object. Later those regions were grouped separately with Ansys CFD-ICEM software.
The microtubule has complex structure and its implementation with mesh generator is not
easy task since each microtubule needs manual generation. Therefore we first tried 2d model
with straight shape of microtubules. This implementation was simply working with spherical
shape of cell model as seen in Figure 16.
(a)
(b)
Figure 16. Simple implementation of straight microtubule structures. Microtubules share the mesh
information with cell meshes. (a) shows the test of active transport by microtubules and (b) is the
1st cell model with consideration of microtubules volume ratio. Microtubules facilitate the vesicle
transport to the cell boundary.
Once we have 3d model from confocal microscopic images, we met the necessary of
curved shape of microtubules mostly spanned from its origin to the end-feet of astrocyte
processes as seen in Figure 1. To imitate more real-like structures without losing the
47
computational compatibility, we developed smooth and shortest path finding algorithm as shown
in section 2.1.2. The algorithm worked to generate smoothly curved microtubules form the
MTOC to the end-feed of each process. All pathways to the end points were kept inside the cell
domain as shown in Figure 17.
Figure 17. Real-like microtubule implementation in the astrocyte 3d cell model. The red curves
represent microtubules.
3.3
Kinetic Model of AQP4 Synthesis
3.3.1 AQP4 Synthesis at Steady State
The steady state of the system describing AQP4 synthesis and transport was matched
with existing data. AQP4 proteins have a known half-life of 24 hours in the cell, providing
information about its degradation kinetics. The abundance of AQP4 mRNA and proteins per cell
48
is assumed to fall within the range of more than four thousand genes measured in mammalian
cells. However, the precise number of transcripts and AQP4 proteins per cell is not known.
Kinetic parameters in the model as shown in eqs. (3) to (8) were determined to match
known data that ~2 copies of aqp4 mRNA and ~900 AQP4 copies per mRNA for an hour at
initial steady state. The shifted transcription rate by SFN treatment was set as ~8 copies of aqp4
mRNA from the increased temporal level of activated DNA. The degradation rate of AQP4 is set
to 24 hour half-life and its corresponding aqp4 mRNA is assumed ~4.8 hours due to the fact that
on average mRNAs are 5 times less stable than the corresponding proteins (Schwanhausser et al.
2011). This information is converted to kinetic rate constants as in Table II.
Table II. Kinetic rates of AQP4 synthesis
Parameters
Values [s-1]
Transcription rate with no SFN treatment
5.55e-4
Transcription rate with the SFN treatment
2.08e-3
Translation rate
0.25
Degradation rate of mRNA
4.01e-5
Degradation rate of AQP4
8.02e-6
3.3.2 Dynamic AQP4 Synthesis
To simulate the dynamics of the system with cellular simulator, we used the proposed
model and explicitly calculated kinetic parameters from previous section. Based on the
observation of AQP4 up-regulation after SFN treatment, our cellular simulator calculates the
relative concentration with the deterministic model and its fluctuation with stochastic simulation.
The simulation result for 18 hours after SFN treatment can be seen in the Figure 18 and
49
conclusively we had AQP4 expression increased by 47.5% after 9 hours and 68.0% after 18
hours. Activated DNA keeps higher level after SFN treatment by induced Nrf2 and returns to
initial level after 8 hours due to the Nrf2 export out from nucleus. We note that mRNA
concentration interestingly decreased after 8 hours by adjusted transcription rate and reduced
Nrf2 for AQP4 gene activation step. A discrete event simulation was performed with the
approximated Gillespie method. This approximation has a discretized time step size of ∆t=1 sec
and the total simulation time is 18 hours. Stochastic simulations were repeated for 100 times and
we confirmed that their average converges to the deterministic simulation as shown in Figure 18.
Even though the stochastic model generates discrete number of molecules, the use of
relative quantities is preferred in the current study since the absolute number of AQP4 per cells is
not known.
Figure 18. Stochastic simulation results of AQP4 synthesis. Solid lines in red color
represent the mean of 100 stochastic simulations and the black dotted lines are the result of
the deterministic model. The proposed model and kinetic parameters match the
experimental results.
50
3.4
Spatial Translocation of Aquaporin-4
In brain astrocytes, AQP4 are highly expressed in the end-feet of the astrocytes. After
folding into vesicles, AQP4 proteins are found to be transported towards the endfeet via the
cytoskeleton. Each vesicle contains numerous AQP4 proteins assembled in arrays. In our
stochastic model, we incorporated two intracellular transport mechanisms; passive transport in
the cytoplasm by pure diffusion and active transport driven by motor protein along microtubules.
Diffusion events of AQP4 proteins to the neighboring mesh volumes take with 0.125 μm 2/s
diffusivity with random direction and convective transport associated to motor protein with the
speed of 1 μm /s in average along microtubules spanned with ~20 μm (Vale et al. 1985; von
Massow et al. 1989; Steinberg and Schliwa 1996; Falk et al. 2009; Falk et al. 2010).
To compare both active and passive trnasport mechanism in terms of distance from
originated point, we tested 2d circular cell models by varing the number of MTs as shown in
Figure 19. In the figure, we can see that the more microtubules in the cell the more it facilitates
Figure 19. Stochastic active and passive transport result with different number of microtubules
during 100 simulation time in the circular 2d cell model. The more microtubules in the cell, the
more it facilitates the active transport to the boundary of cell model which is assumed plasma
membrane in the real 3d cell model. In the simulation of 1000 molecules, there are 21 molecules
for the case of the number of 5 MTs, 93 molecules for 100 MTs and ~390 molecules for 300 MTs
associated with active transport by MTs.
51
the active transport to the boundary of cell model which is assumed plasma membrane in the real
3d cell model. In the simulation of 1000 molecules, there are 21 molecules for the case of the
number of 5 MTs, 93 molecules for 100 MTs and ~390 molecules for 300 MTs associated with
active transport by MTs. In terms of distance, this model has a radius, 100 μm, from the center to
boundary and average distances of entire molecules of all three cases are respectively 11.7 μm,
15.5 μm, and 28.7 μm. This comparison confirmed that only pure diffusion does not transport the
molecules to the cell membrane as well as transportation to the targeted region is facilitated by
active transport along microtubules. But due to the extremely simplified cell domain and
microtubules with shape of straight line, the difference of traveling distance is only effective as
in this comparison.
To accomplish the real-like spatial intracellular transport of AQP4, we developed 3d cell
model from confocal microscopies which captured curtured astrocyte. In that model, the
microtubule network originating from the centrosome is built by an efficient ahort and smooth
path-finding algorithm developed by the author and the stochastic simulation result of AQP4
spatiotermporal intracellular transport in the astrocyte is shown in Figure 20. In Figure 20 (b),
increased amount of AQP4 (about 68% after 18 hours) is mostly transported to the end-feet of
processes. Each dot in Figure 20 (b) has different number of AQP4 molecules from 1 to 50 by
random formation and about 800 vesicles (corresponding to ~2×105 AQP4 molecules) were
formed from ER membrane. In Figure 20 (a), green dots stand for AQP4 vesicles after 18 hour
and those are mostly concentrated at the end-feet of each process but there are still found a
number of dots throughout the cytoplasm which is now on the microtubule tracks.
52
10μm
(a)
(b)
Figure 20. AQP4 in primary cultured astrocytes (a) and model (b). In (a), AQP4 vesicles are
stained as green dot and blue circle in the center is nucleus. AQP4 are seen in vesicles as
well as along astrocytic processes during their transport towards the endfeet. Scale bar in
both images indicates 10 microns. The computational model of astrocyte in (b) visualizes
realistic microtubules and size-varied blue dot AQP4 vesicles in 3d distributed system.
3.5
Nrf2 activation by SFN treatment in astrocytes
In this section, we summarized the Nrf2 localization experiment result. The up-regulated
localization of Nrf2 was used in the cellular simulator with careful interpretation. For quantifying
the nuclear translocation of phosphorylated Nrf2, cells were treated with SFN for 15min, 30min,
1.5h, 3h, 9h, and 18h and the fluorescence intensity inside the nuclear region was normalized
with that of DMSO-vehicle control. The purity of the astrocyte culture was confirmed using
astrocyte-specific marker glial fibrillary acidic protein (GFAP), image not shown. The timedependent increase and eventual normalization of pNrf2 levels inside the nucleus is shown in
Figure 21.
53
Figure 21. Translocation of Nrf2 from cytoplasm to nucleus after SFN exposure. In
control cells, Nrf2 is found in cytoplasm with a weak basal immunoreactivity in the
nucleus. After SFN exposure, Nrf2 is translocated to the nucleus showing a high
immunoreactivity, and Nrf2-filled vesicles marking transport of Nrf2 are observed.
Nuclear intensity of pNrf2 is increased by more than 60% by the 15min time point, in
agreement with findings by Jain et al. stating that Nrf2 translocation occurs as early as 15
minutes (Jain et al. 2005). The initial increase in pNRf2 nuclear intensity is followed by a
decrease, likely due to the depletion of the pNrf2 pool in the nucleus after de-phosphorylation
and nuclear export of Nrf2. pNrf2 levels show normalization by the 9 hour time point in Figure
21. The time point of pNrf2 normalization after treatment is in agreement with previous studies.
A study by Jain et al. reports that NRf2 starts to exit the nucleus between 1-4 hours after
treatment and achieves normal levels at 8 hours (Jain et al. 2005).
3.6
Up-regulation of AQP4 by SFN
Western blot (WB) was performed to quantify the up-regulation of AQP4 by 13µM of
SFN. For SDS-PAGE gel electrophoresis, cell lysates were prepared after 9 and 18 hours of
continuous SFN exposure, and Bradford assay was performed to ensure equal loading of proteins
54
in all samples. The western blot data was analyzed with densitometry in Image J (Mao et al.
2006). The results show that AQP4 expression was 1.48 and 1.68 fold higher after 9 and 18
hours of SFN exposure as shown in Figure 22. The data confirmed that SFN induces upregulation of AQP4 in astrocyte cells, in agreement with previous findings that SFN injection
induces AQP4 up-regulation in the brain. Zhao et al find a 65% AQP4 up-regulation by SFN
administration for 24 hours to be effective in an in vivo model of traumatic brain injury (Zhao et
al. 2005).
The modeling of cellular processes leading to AQP4 up-regulation will be presented in
the next section. The cellular simulator can identify kinetic parameters that agree with
experimentally measured trends.
Figure 22. Western blotting shows that SFN induced a 48% up-regulation after 9 hours of
continuous exposure, and 68% up-regulation after 18 hours compared to control.
55
4
DISCUSSION
The present work aims at developing a detailed mathematical modeling and simulation of
AQP4 intracellular transport, which allows modeling of the spatiotemporal dynamics of cellular
level transport in 3d and coupling of reaction and diffusion process with stochastic fluctuation
involved in each reaction and diffusion procedure. As long as we simply the detailed complex
mechanism with lumped kinetic model and approximate the stochastic dynamics avoiding the
record of all possible scenarios in terms with probability distribution for a single simulation, our
approach is still feasible. In addition, the framework of reaction and diffusion stochastic
simulation is based on multi-particle based simulation benefits the successful outcomes.
Although there is currently only single case study in this thesis, we believe this approach will be
applied to interpret other proteins’ intracellular transport mechanism as well.
Klann et al. proposed the first vesicle transport model tracking individual (agent based)
vesicles through the complete cell (Klann et al. 2009). While their work also used approximated
stochastic simulation in the spherical shape of cell model and unrealistic cytoskeleton structure,
the sequential actions of vesicle transport nevertheless satisfied faithfully principles of vesicle
transport true functions. In our present work, we improved the cell model with true model from
cultured astrocyte even though its intracellular environment is simplified with a couple of
compartment units which are necessary for AQP4 transport mechanism by lumping all other
kinetics. Moreover, we are first to result in both quantitative and qualitative identification of
vesicle transport in the cellular level with computational approach. Especially, the
56
implementation of real-like microtubule structure increased the value of the originality of current
work.
To support our model and provide unknown parameter of Nrf2 concentration, we had
experiments of AQP4 quantification in the temporal response of a cell towards a disturbance or a
therapeutic molecule. This analysis might facilitate our understanding of SFN dose-response to
Nrf2 and AQP4 on an organism-level, the present mathematical model in this thesis takes most
of kinetic parameters from existing transcription and translation kinetics. Moreover the upregulated simulation result after 9 and 18 hours matched with existing AQP4 up-regulation
experimental results is meaningful accomplishment. Moreover the approximated Gillespie
algorithm reduced the computational cost and approach reasonable expectation with stochasticity
by calculating sum of events during discrete time intervals. Poisson distribution used to the
sampling of random numbers from the probability during reasonably short period ∆t properly
applied to give positive number of particle involved in the reaction mechanism. However,
therapeutically effective SFN dosage to achieve significant AQP4 up-regulation is beyond of the
scope of this study and the determination of in-silico AQP4 protein up-regulation should be
addressed with full understanding on the mechanism of AQP4 synthesis and spatiotemporal
transport and tons of experimental efforts.
So far, the computational effort of particle based simulations and reasonable
approximation in the mathematical model does allow the handling of real-like micro-scale
dynamics in a single cell. The present simulation is also pioneering steps to understand
biological system with computational solution and 3d visualization. However, it is true that our
model is still far from the real intracellular space and the up-regulation outcome with Nrf2
57
increases by SFN administration might be obtained superficially. This model may or may not
address useful pathological control strategy corresponding to some type of cerebral edema in
terms of water homeostasis, but we cannot but agree that our approach has too many things to be
improved in aspects of validity, reality, clarity, and etc.
58
5
CONCLUSION
Systems biology is now one of emerging trends in the biology study. Its quantitative
modeling assisted by mathematical model and computational simulation on the neurophysiology
is now capable to represent the microscale mechanisms. If those quantitative information about
the dynamics of protein synthesis can be combined with the existing interpretation on the system,
it must be immensely useful information about the life of a protein in a cell, from its generation,
transport and expression to its degradation.
In this study, the entire 'life-story' of a protein including the transcriptional activation of
its gene, the transcription of corresponding mRNA and transport of protein from the nucleus, and
its trafficking to the destination and eventual recycling was re-addressed from mathematical
model. The dynamics of protein trafficking and transport are modeled and trajectories of proteins
in a cell level were visualized by solving the proposed mathematical model. Even though understanding biological systems such as metabolic networks or cell signaling networks still has a
massive limitation, still our effort to integrate computational approach and biological systems is
highly advantageous and its further application must be supported by this work.
59
CITED LITERATURE
Agre, P., S. Sasaki, et al. (1993). "Aquaporins: a family of water channel proteins." American
Journal of Physiology - Renal Physiology 265(3): F461.
Alikina, T. Y., N. B. Illarionova, et al. (2012). "Identification of new M23A mRNA of mouse
aquaporin-4 expressed in brain, liver, and kidney." Biochemistry (Moscow) 77(5): 425434.
Ansys (2013). Ansys ICEM CFD.
Arkin, A., J. Ross, et al. (1998). "Stochastic Kinetic Analysis of Developmental Pathway
Bifurcation in Phage λ-Infected Escherichia coli Cells." Genetics 149(4): 1633-1648.
Auger, A., P. Chatelain, et al. (2006). "R-leaping: Accelerating the stochastic simulation
algorithm by reaction leaps." The Journal of Chemical Physics 125(8): 084103.
Badaut, J., F. Lasbennes, et al. (2002). "Aquaporins in Brain[colon] Distribution, Physiology,
and Pathophysiology." J Cereb Blood Flow Metab 22(4): 367-378.
Baras, F. and M. M. Mansour (1996). "Reaction-diffusion master equation: A comparison with
microscopic simulations." Physical Review E 54(6): 6139-6148.
Beisson, J. and M. Wright (2003). "Basal body/centriole assembly and continuity." Current
Opinion in Cell Biology 15(1): 96-104.
Benga, G., O. Popescu, et al. (1986). "Water permeability in human erythrocytes: identification
of membrane proteins involved in water transport." European journal of cell biology
41(2): 252-262.
Bloch, O., K. I. Auguste, et al. (2006). "Accelerated progression of kaolin-induced
hydrocephalus in aquaporin-4-deficient mice." J Cereb Blood Flow Metab 26(12): 15271537.
Böhm, K. J., P. Steinmetzer, et al. (1997). "Kinesin-driven microtubule motility in the presence
of alkaline–earth metal ions: Indication for a calcium ion-dependent motility." Cell
Motility and the Cytoskeleton 37(3): 226-231.
Bruce Alberts, A. J., Julian Lewis, Martin Raff, Keith Roberts, Peter Walter (2002). Molecular
Biology of the Cell, New York: Garland Science.
Cai, X. and Z. Xu (2007). "K-leap method for accelerating stochastic simulation of coupled
chemical reactions." The Journal of Chemical Physics 126(7): 074102.
Cho, H.-Y., S. P. Reddy, et al. (2005). "Gene expression profiling of NRF2-mediated protection
against oxidative injury." Free Radical Biology and Medicine 38(3): 325-343.
Crane, J. M., J. L. Bennett, et al. (2009). "Live Cell Analysis of Aquaporin-4 M1/M23
Interactions and Regulated Orthogonal Array Assembly in Glial Cells." Journal of
Biological Chemistry 284(51): 35850-35860.
Darzacq, X., Y. Shav-Tal, et al. (2007). "In vivo dynamics of RNA polymerase II transcription."
Nat Struct Mol Biol 14(9): 796-806.
Denker, B. M., B. L. Smith, et al. (1988). "Identification, purification, and partial
characterization of a novel Mr 28,000 integral membrane protein from erythrocytes and
renal tubules." Journal of Biological Chemistry 263(30): 15634-15642.
Ellis, R. J. (2001). "Macromolecular crowding: an important but neglected aspect of the
intracellular environment." Current Opinion in Structural Biology 11(1): 114-119.
60
Elowitz, M. B., M. G. Surette, et al. (1999). "Protein mobility in the cytoplasm of Escherichia
coli." J Bacteriol 181(1): 197-203.
Feldman, J. L., S. Geimer, et al. (2007). "The Mother Centriole Plays an Instructive Role in
Defining Cell Geometry." PLoS Biol 5(6): e149.
Fenton, R., H. Moeller, et al. (2010). "Differential water permeability and regulation of three
aquaporin-4 isoforms." Cellular and Molecular Life Sciences 67(5): 829-840.
Finnie, J. W., J. Manavis, et al. (2008). "Aquaporin-4 in Acute Cerebral Edema Produced by
Clostridium perfringens Type D Epsilon Toxin." Veterinary Pathology Online 45(3):
307-309.
Frigeri, A., M. A. Gropper, et al. (1995). "Localization of MIWC and GLIP water channel
homologs in neuromuscular, epithelial and glandular tissues." Journal of Cell Science
108(9): 2993-3002.
Gibson, M. A. and J. Bruck (2000). "Efficient Exact Stochastic Simulation of Chemical Systems
with Many Species and Many Channels." The Journal of Physical Chemistry A 104(9):
1876-1889.
Gillespie, D. T. (1977). "Exact stochastic simulation of coupled chemical reactions." The Journal
of Physical Chemistry 81(25): 2340-2361.
Gillespie, D. T. (2001). "Approximate accelerated stochastic simulation of chemically reacting
systems." The Journal of Chemical Physics 115(4): 1716-1733.
Hasegawa, H., T. Ma, et al. (1994). "Molecular cloning of a mercurial-insensitive water channel
expressed in selected water-transporting tissues." Journal of Biological Chemistry 269(8):
5497-5500.
Hazama, A., D. Kozono, et al. (2002). "Ion Permeation of AQP6 Water Channel Protein:
SINGLE-CHANNEL RECORDINGS AFTER Hg2+ACTIVATION." Journal of
Biological Chemistry 277(32): 29224-29230.
Hirt, L., B. Ternon, et al. (2009). "Protective role of early Aquaporin-4 induction against
postischemic edema formation." Journal of Cerebral Blood Flow and Metabolism 29(2):
423-433.
Hong, W. (1998). "Protein transport from the endoplasmic reticulum to the Golgi apparatus." J.
Cell Sci. 111: 2831-2839.
Isaacson, S. A. and D. Isaacson (2009). "Reaction-diffusion master equation, diffusion-limited
reactions, and singular potentials." Phys Rev E Stat Nonlin Soft Matter Phys 80(6 Pt 2):
066106.
Jain, A. K., D. A. Bloom, et al. (2005). "Nuclear Import and Export Signals in Control of Nrf2."
Journal of Biological Chemistry 280(32): 29158-29168.
Jin, B.-J., A. Rossi, et al. (2011). "Model of Aquaporin-4 Supramolecular Assembly in
Orthogonal Arrays Based on Heterotetrameric Association of M1-M23 Isoforms."
Biophysical Journal 100(12): 2936-2945.
Jung, J. S., R. V. Bhat, et al. (1994). "Molecular characterization of an aquaporin cDNA from
brain: candidate osmoreceptor and regulator of water balance." Proceedings of the
National Academy of Sciences 91(26): 13052-13056.
Kang, M. and H. G. Othmer (2009). "Spatiotemporal characteristics of calcium dynamics in
astrocytes." Chaos: An Interdisciplinary Journal of Nonlinear Science 19(3): 037116037121.
61
Ke, C., W. S. Poon, et al. (2001). "Heterogeneous responses of aquaporin-4 in oedema formation
in a replicated severe traumatic brain injury model in rats." Neuroscience Letters 301(1):
21-24.
Kerr, I. M., A. P. Costa-Pereira, et al. (2003). "Of JAKs, STATs, blind watchmakers, jeeps and
trains." FEBS Lett 546(1): 1-5.
Khanin, R. and D. J. Higham (2008). "Chemical Master Equation and Langevin regimes for a
gene transcription model." Theoretical Computer Science 408(1): 31-40.
Kholodenko, B. N. (2003). "Four-dimensional organization of protein kinase signaling cascades:
the roles of diffusion, endocytosis and molecular motors." J Exp Biol 206(Pt 12): 20732082.
Kiening, K. L., F. K. H. van Landeghem, et al. (2002). "Decreased hemispheric Aquaporin-4 is
linked to evolving brain edema following controlled cortical impact injury in rats."
Neuroscience Letters 324(2): 105-108.
Kimelberg, H. K. (2004). "Water homeostasis in the brain: Basic concepts." Neuroscience 129(4):
851-860.
Kobayashi, A., M.-I. Kang, et al. (2004). "Oxidative Stress Sensor Keap1 Functions as an
Adaptor for Cul3-Based E3 Ligase To Regulate Proteasomal Degradation of Nrf2."
Molecular and Cellular Biology 24(16): 7130-7139.
Lee, P. S., L. B. Shaw, et al. (2003). "Insights into the relation between mrna and protein
expression patterns: ii. Experimental observations in Escherichia coli1." Biotechnology
and Bioengineering 84(7): 834-841.
Li, B., S. Chen, et al. (2012). "Modeling the Contributions of Ca<sup>2+</sup> Flows to
Spontaneous Ca<sup>2+</sup> Oscillations and Cortical Spreading DepressionTriggered Ca<sup>2+</sup> Waves in Astrocyte Networks." PLoS ONE 7(10): e48534.
Lu, M., M. D. Lee, et al. (1996). "The human AQP4 gene: definition of the locus encoding two
water channel polypeptides in brain." Proceedings of the National Academy of Sciences
93(20): 10908-10912.
Luby-Phelps, K. (1999). Cytoarchitecture and Physical Properties of Cytoplasm: Volume,
Viscosity, Diffusion, Intracellular Surface Area. International Review of Cytology. D. E.
B. Harry Walter and A. S. Paul, Academic Press. 192: 189-221.
MacDonald, C. T., J. H. Gibbs, et al. (1968). "Kinetics of biopolymerization on nucleic acid
templates." Biopolymers 6(1): 1-25.
Mao, X., T. L. Enno, et al. (2006). "Aquaporin-4 changes in rat brain with severe
hydrocephalus." European Journal of Neuroscience 23(11): 2929-2936.
Materialise Mimics.
MathWorks (2013). MATLAB. R2012b.
McAllister, J. P. and J. M. Miller (2006). "Aquaporin-4 and hydrocephalus." Journal of
Neurosurgery: Pediatrics 105(6): 457-458.
McMahon, M., K. Itoh, et al. (2003). "Keap1-dependent Proteasomal Degradation of
Transcription Factor Nrf2 Contributes to the Negative Regulation of Antioxidant
Response Element-driven Gene Expression." Journal of Biological Chemistry 278(24):
21592-21600.
Mehra, A. and V. Hatzimanikatis (2006). "An Algorithmic Framework for Genome-Wide
Modeling and Analysis of Translation Networks." Biophysical Journal 90(4): 1136-1146.
62
Moe, S. E., J. G. Sorbo, et al. (2008). "New isoforms of rat Aquaporin-4." Genomics 91(4): 367377.
Morishita, Y., T. Matsuzaki, et al. (2005). "Disruption of Aquaporin-11 Produces Polycystic
Kidneys following Vacuolization of the Proximal Tubule." Molecular and Cellular
Biology 25(17): 7770-7779.
Müller, M. J. I., S. Klumpp, et al. (2008). "Tug-of-war as a cooperative mechanism for
bidirectional cargo transport by molecular motors." Proceedings of the National Academy
of Sciences 105(12): 4609-4614.
Nagelhus, E. A., M. L. Veruki, et al. (1998). "Aquaporin-4 Water Channel Protein in the Rat
Retina and Optic Nerve: Polarized Expression in Müller Cells and Fibrous Astrocytes."
The Journal of Neuroscience 18(7): 2506-2519.
Nedergaard, M., B. Ransom, et al. (2003). "New roles for astrocytes: Redefining the functional
architecture of the brain." Trends in Neurosciences 26(10): 523-530.
Nielsen, S., E. Arnulf Nagelhus, et al. (1997). "Specialized Membrane Domains for Water
Transport in Glial Cells: High-Resolution Immunogold Cytochemistry of Aquaporin-4 in
Rat Brain." The Journal of Neuroscience 17(1): 171-180.
Nuriya, M., T. Shinotsuka, et al. (2012). "Diffusion Properties of Molecules at the Blood–Brain
Interface: Potential Contributions of Astrocyte Endfeet to Diffusion Barrier Functions."
Cerebral Cortex.
Paul, L., M. Madan, et al. (2009). "The altered expression of aquaporin 1 and 4 in choroid plexus
of congenital hydrocephalus." Cerebrospinal Fluid Research 6(Suppl 1): 1-1.
Paulsson, J., O. G. Berg, et al. (2000). "Stochastic focusing: fluctuation-enhanced sensitivity of
intracellular regulation." Proc Natl Acad Sci U S A 97(13): 7148-7153.
Preston, G. M., T. P. Carroll, et al. (1992). "Appearance of Water Channels in Xenopus Oocytes
Expressing Red Cell CHIP28 Protein." Science 256(5055): 385-387.
Rodríguez, J. V., J. A. Kaandorp, et al. (2006). "Spatial stochastic modelling of the
phosphoenolpyruvate-dependent phosphotransferase (PTS) pathway in Escherichia coli."
Bioinformatics 22(15): 1895-1901.
Ross, J. L., M. Y. Ali, et al. (2008). "Cargo transport: molecular motors navigate a complex
cytoskeleton." Current Opinion in Cell Biology 20(1): 41-47.
Rossi, A., J. M. Crane, et al. (2011). "Aquaporin-4 Mz isoform: Brain expression,
supramolecular assembly and neuromyelitis optica antibody binding." Glia 59(7): 10561063.
Saadoun, S., M. C. Papadopoulos, et al. (2002). "Aquaporin-4 expression is increased in
oedematous human brain tumours." Journal of Neurology, Neurosurgery & Psychiatry
72(2): 262-265.
Schwanhausser, B., D. Busse, et al. (2011). "Global quantification of mammalian gene
expression control." Nature 473(7347): 337-342.
Seifert, G., K. Schilling, et al. (2006). "Astrocyte dysfunction in neurological disorders: a
molecular perspective." Nat Rev Neurosci 7(3): 194-206.
Shen, X. Q., M. Miyajima, et al. (2006). "Expression of the water-channel protein aquaporin-4 in
the H-Tx rat: possible compensatory role in spontaneously arrested hydrocephalus."
Journal of Neurosurgery: Pediatrics 105(6): 459-464.
63
Silberstein, C., R. Bouley, et al. (2004). "Membrane organization and function of M1 and M23
isoforms of aquaporin-4 in epithelial cells." American Journal of Physiology - Renal
Physiology 287(3): F501-F511.
Simard, M. and M. Nedergaard (2004). "The neurobiology of glia in the context of water and ion
homeostasis." Neuroscience 129(4): 877-896.
Skjolding, A., I. Rowland, et al. (2010). "Hydrocephalus induces dynamic spatiotemporal
regulation of aquaporin-4 expression in the rat brain." Cerebrospinal Fluid Research 7(1):
20.
Sofroniew, M. and H. Vinters (2010). "Astrocytes: biology and pathology." Acta
Neuropathologica 119(1): 7-35.
Solenov, E., H. Watanabe, et al. (2004). "Sevenfold-reduced osmotic water permeability in
primary astrocyte cultures from AQP-4-deficient mice, measured by a fluorescence
quenching method." American Journal of Physiology - Cell Physiology 286(2): C426C432.
Steinberg, G. and M. Schliwa (1996). "Characterization of the Biophysical and Motility
Properties of Kinesin from the Fungus Neurospora crassa." Journal of Biological
Chemistry 271(13): 7516-7521.
Strand, L., S. E. Moe, et al. (2009). "Roles of Aquaporin-4 Isoforms and Amino Acids in Square
Array Assembly." Biochemistry 48(25): 5785-5793.
Suzuki, R., M. Okuda, et al. (2006). Astrocytes co-express aquaporin-1, -4, and vascular
endothelial growth factor in brain edema tissue associated with brain contusion. Brain
Edema XIII. J. Hoff, R. Keep, G. Xi and Y. Hua, Springer Vienna. 96: 398-401.
Tajima, M., J. M. Crane, et al. (2010). "Aquaporin-4 (AQP4) Associations and Array Dynamics
Probed by Photobleaching and Single-molecule Analysis of Green Fluorescent ProteinAQP4 Chimeras." Journal of Biological Chemistry 285(11): 8163-8170.
Tomás-Camardiel, M., J. L. Venero, et al. (2004). "In vivo expression of aquaporin-4 by reactive
microglia." Journal of Neurochemistry 91(4): 891-899.
Tomás, M., P. Marín, et al. (2005). "Ethanol perturbs the secretory pathway in astrocytes."
Neurobiology of Disease 20(3): 773-784.
Tourdias, T., I. Dragonu, et al. (2009). "Aquaporin-4 correlates with apparent diffusion
coefficient and hydrocephalus severity in the rat brain: A combined MRI–histological
study." NeuroImage 47(2): 659-666.
Trcek, T., Daniel R. Larson, et al. (2011). "Single-Molecule mRNA Decay Measurements
Reveal Promoter- Regulated mRNA Stability in Yeast." Cell 147(7): 1484-1497.
Trinh, S., P. Arce, et al. (2000). "Effective Diffusivities of Point-Like Molecules in Isotropic
Porous Media by Monte Carlo Simulation." Transport in Porous Media 38(3): 241-259.
Umenishi, F. and A. S. Verkman (1998). "Isolation and Functional Analysis of Alternative
Promoters in the Human Aquaporin-4 Water Channel Gene." Genomics 50(3): 373-377.
Vale, R. D., T. S. Reese, et al. (1985). "Identification of a novel force-generating protein, kinesin,
involved in microtubule-based motility." Cell 42(1): 39-50.
Van Kampen, N. G. (2011). Stochastic Processes in Physics and Chemistry, Elsevier Science.
Venero, J. L., M. a. L. Vizuete, et al. (2001). "Aquaporins in the central nervous system."
Progress in Neurobiology 63(3): 321-336.
Verkman, A. S. (2011). "Aquaporins at a glance." Journal of Cell Science 124(13): 2107-2112.
64
Verkman, A. S., D. K. Binder, et al. (2006). "Three distinct roles of aquaporin-4 in brain function
revealed by knockout mice." Biochimica et Biophysica Acta (BBA) - Biomembranes
1758(8): 1085-1093.
Voeltz, G. K., M. M. Rolls, et al. (2002). "Structural organization of the endoplasmic reticulum."
EMBO Rep 3(10): 944-950.
von Massow, A., E. M. Mandelkow, et al. (1989). "Interaction between kinesin, microtubules,
and microtubule-associated protein 2." Cell Motil Cytoskeleton 14(4): 562-571.
Vossenkämper, A., P. I. Nedvetsky, et al. (2007). "Microtubules are needed for the perinuclear
positioning of aquaporin-2 after its endocytic retrieval in renal principal cells." American
Journal of Physiology - Cell Physiology 293(3): C1129-C1138.
Wang, F., N. A. Smith, et al. (2012). "Astrocytes Modulate Neural Network Activity by Ca2+Dependent Uptake of Extracellular K+." Sci. Signal. 5(218): ra26-.
Zador, Z., S. Stiver, et al. (2009). Role of Aquaporin-4 in Cerebral Edema and Stroke.
Aquaporins. E. Beitz, Springer Berlin Heidelberg. 190: 159-170.
Zelenin, S., E. Gunnarson, et al. (2000). "Identification of a New Form of AQP4 mRNA That Is
Developmentally Expressed in Mouse Brain." Pediatr Res 48(3): 335-339.
Zhao, J., A. N. Moore, et al. (2005). "Sulforaphane enhances aquaporin-4 expression and
decreases cerebral edema following traumatic brain injury." Journal of Neuroscience
Research 82(4): 499-506.
PART II. DESIGN AND MODELING OF
COMPLEX DISTILLATION PROCESS
ENERGY
EFFICIENT
Of the many separation technologies in use, distillation continues to be the most widely
used and versatile method of separation, particularly in the petrochemical industry. Separation
processes make up 40%-70% of capital and operating costs of chemical manufacturing.
Accordingly, even a small improvement in the distillative separations would yield huge energy
savings. To harness possible energy improvements, the use of thermodynamically integrated
complex columns has been suggested (Doherty and Malone 2001; Engelien and Skogestad 2005;
Grossmann et al. 2005; Taylor 2007; Harwardt et al. 2008; Holland et al. 2009; Giridhar and
Agrawal 2010). Moreover, the rising cost of energy and increased concern over atmospheric
carbon emission has sparked interest in the further optimization of distillation systems. Recent
work by Luyben’s group (Ling and Luyben 2009) shows practical and effective control strategies
for heat integrated columns including divided wall column which was first proposed by Wright
(Wright 1949).
Therefore, we aim at novel methods for discovering energy-efficient configurations
systematically in this part. Three or more columns are often necessary to sufficiently separate
components of a mixture and different configurations of columns in a network can achieve the
same separation with smaller amount of energy. Until recently, however, there were no rigorous
algorithms to synthesize complete separation flowsheets without limiting assumptions of the
thermodynamic vapor-liquid equilibrium model or restrictions in its applicability like sharp splits
or nondistributing species. Albeit some structural approaches specifically addressing
synthesizing complex network structures have led to the possibility of enumerating all possible
configurations systematically (Agrawal 2003; Giridhar and Agrawal 2010), a detailed design
65
66
including realizable column profiles is needed to judge whether a certain configuration can
actually achieve the desired purity.
We extended the minimum bubble point distance algorithm and temperature collocation
methodology to rigorously design complex networks to separate ideal and nonideal
multicomponent mixtures into products of desired purity using heat-integrated prefractionating
columns. We also employed inverse design procedure enables the systematic design of
physically realizable separations for mixtures with a large number of species. The computer
procedure robustly converges to the desired purity targets, unless the desired purity target is
thermodynamically impossible to realize. The algorithm also rapidly identifies infeasible
specifications without fail. Finally, synthesized networks were validated with AspenPlus
matching exactly the inverse design results with the target purity. The rigorous flowsheet design
combined with validation of the networks with commercial flowsheet simulators enables the
systematic design of energy-efficient separation networks. The methodology is ready to address
currently unresolved design problems such as the computer-aided design of energy-efficient
separations, the design of biorefineries, or new process designs for carbon sequestration.
67
1
INTRODUCTION TO THE DESIGN AND SYNTHESIS OF COMPLEX
DISTILLATION PROCESS
Distillation continues to be the most widely used separation technique in the chemical
and petrochemical industry. Although distillation requires heat, it still remains a simple and
economic way to separate many multicomponent liquid mixtures. Despite the common
misconception of distillation as a mature field, complex column configurations for enhanced
energy-efficiency, extractive, or reactive separations, are an open field of research (Siirola 1996;
Widagdo and Seider 1996; Malone and Doherty 2000; Barnicki and Siirola 2004). Recently,
advancements in design and synthesis of separation systems slowly do away with this myth.
Doherty and coworkers (Van Dongen and Doherty 1985; Doherty and Malone 2001) developed
the boundary value method based on the intersection of stripping and rectifying profiles in the
continuous stage number space. Marquardt(Marquardt et al. 2008) proposed the use of rectifying
bodies, which are constructed by connecting the pinch points belonging to column sections.
According to this method, intersection of the rectifying bodies between two adjacent rectifying
and stripping column sections implies feasibility. Yet, realizable column profiles are not
necessarily confined in the convex hull spanned by rectifying bodies, nor do pinch points always
come to lie in real space so as to form a convex geometric object (Holland et al. 2004).
In the synthesis and design of separation networks, Agrawal and Fidkowski(Agrawal and
Fidkowski 1998) proposed a strategy to develop structural alternatives for complex and simple
distillation networks. Another important advancement is Agrawal’s procedure for the exhaustive
enumeration of complex flowsheets (Agrawal 2003). They systematically evaluated the energy
68
consumption of so called basic feasible configurations of complex networks and advocate the
energy savings that can potentially be achieved through feasible thermal coupling. Their
complete synthesis of candidate flowsheet structures is an important element in a fully automatic
separation synthesis methodology.
Lucia and Taylor (Lucia and Taylor 2006; Lucia and Taylor 2007) studied composition
trajectories related to stationary points and residual composition maps to understand the
separation path in different boundary regions of composition planes. Recently, they propose a
new approach based on the concept of shortest stripping line distance to identify minimum
energy requirements by minimizing the length of the stripping profile (Lucia et al. 2008).
Glasser and Hildebrandt (Holland et al. 2004; Tapp et al. 2004) developed a novel design
technique using generalizing column section computations giving rise to column profile maps
encompassing all realizable liquid composition profiles for simple as well as complex column
sections. This concept is closely related to residue curve maps and can be applied for assessing
feasible separations.
Our group developed a novel composition model approach called temperature collocation,
replacing the conventional tray number by the bubble point temperature (Zhang and Linninger
2004). This thermodynamic variable transformation drastically reduces the design search space,
because this absolute bubble point temperature is independent of stage numbers and applies to
both stripping and rectifying sections. Every new stage computation which traditionally opens
new search dimensions for its number of stages or section length can be transformed onto a
global and characteristic temperature scale. This efficient and reliable computational algorithm
entitled Bubble Point Distance criterion, BPD, can establish the feasibility or impossibility of
69
realizing desired product specifications for sharp or sloppy separations of multicomponent
mixtures using ideal, nonideal or azeotropic vapor-liquid equilibrium models. The BPD criterion
was previously applied successfully for the automatic synthesis of simple column networks
(Zhang and Linninger 2006; Ruiz et al. 2009).
70
2
METHODOLOGIES OF THE COMPUTATIONAL MODELING OF COMPLEX
COLUMN CONFIGURATION
2.1
Vapor-Liquid Phase Equilibrium Model
In all stages of the column there is a vapor and liquid contact, which generates a mass
transfer phenomenon. Although several vapor-liquid equilibrium models appear to describe the
phase relationship in the many literatures, there is no universal model in nonideal mixtures.
According to the systems and purposes, however, particular VLE models will be better one
rather than others. In this section we will illustrate two general models useful for ideal mixtures;
(1) the constant volatility model and (2) Raoult’s law model. Both models are often close except
for mixtures having very close or very different boiling points. For nonideal mixtures both
models can be applied for the special limiting cases.
2.1.1 Constant Relative Volatility Model
Until recently, the constant relative volatility model is widely used in chemical
engineering distillation processes. The relative volatility (α) of the component is measure of
vapor pressure of the component which also indicates the ease or difficulty to separate in a
mixture. Therefore it is assumed as constant independent of changes in composition and
temperature. Equations (26) and (27) describe the vapor-liquid equilibrium (VLE) in multicomponents mixture.
yi 
 i xi
n

j 1
j
xj
(26)
71
xi 
yi  i
n
y
j 1
j
(27)
j
where i and j are the indices of components, αi (or αj) is the relative volatility of component i
(or j), x is the liquid phase molar fraction, and y is the vapor phase molar fraction. Relative
volatility model is useful for the separation process involved with vapor-liquid equilibrium
stages but may not be applied to the separation in which some components react with each other.
2.1.2 Raoult’s Law Model
Raoult’s law is another model to describe the thermodynamic vapor-liquid equilibrium. It
uses the boiling point of the pure chemical substances. To be applicable, the mixture must be
assumed that the vapor phase is an ideal gas and the liquid phase is an ideal solution. Equation
(28) represents the VLE with Raoult’s law;
yi P  xi Pi sat
where P is the total pressure of the system and
(28)
is the vapor pressure of pure species i at the
temperature T of the system. The Antoine equation is the more generally used expression to
calculate
,
log  Pi sat   A 
B
T C
(29)
where A, B, and C are constants available for a large number of species. Without assuming the
liquid phase is an ideal solution which considers the activity coefficient (γi) into eq. (28), the
Raoult’s law becomes a much more realistic equation for VLE:
yi P  xi i Pi sat
(30)
72
2.2
Complex Column Model
Theoretical work under the minimum vapor flow criteria has shown that complex column
configurations can save up to 70% more energy than conventional column arrays (Hilde K.
Engelien 2005). A typical complex column is the arrangement shown in Figure 23. This
configuration has two feed stream and three product streams; it is divided in four sections: two
rectifying and two stripping profiles. A design is feasible if and only if the composition profiles
between the neighboring two sections intersect.
A
AB
B
ABC
C
BC
(a)
(b)
Figure 23. A complex column network with two feed streams and three product streams. (a)
depicts skeptical configuration of simple and complex column configuration equivalent to
Petlyuk column and (b) shows the sectional decomposition of complex column with
implemented nomenclature of vapor (Vi), liquid (Li), concentration (xi), and Feed (Fj) for each
section i or feed input j.
73
2.2.1 Component Balances for Complex Columns
To get a fully specified system and compute the composition profiles of all complex
sections, first we have to solve a linear algebraic equation system of mass balances. By the same
way we are using intensive and extensive properties to raise and solve the mass balance
equations.
Intensive Properties
n1F f1FA  n2 F f 2 FA  n3F f 3FA  FA
(31)
n1F f1FB  n2 F f 2 FB  n3F f 3FB  FB
(32)
n1F f1FC  n2 F f 2 FC  n3F f 3FC  FC
(33)
f1FA  f1FB  f1FC  1
(34)
f 2 FA  f 2 FB  f 2 FC  1
(35)
f 3FA  f 3FB  f 3FC  1
(36)
s  1  q FA

rq
FC
(37)
We have seven equations eq. (31) to eq. (37) in function of the net component feed flows
2
( n1F , n2 F , n3 F , where ni F   xi Fk Fk
k 1
where xi
Fj
i  1, 2, 3
), the fraction recovery defined as fi j  Fj xi
F
Fj
ni F
is the molar composition of component i at the product stream j ( j  A, B, C ), the reflux
ratio ( r ), and the feed quality ( q ). Then variables – equation: 7 – 7 = 0, now the system is
completely specified. This system of equations can be solved easily expressing eqs. (31) to (37)
as vector and matrix arrangements:
74
 f1FC  1  f1FA   f1FB 
 FC     FA   FB 
 f 2   1   f 2    f 2 
 f 3FC  1  f 3FA   f 3FB 

 



 FA   f1 A
 F    f FB
 B  1
 FC   f1FC
F
f 2 FA
f 2 FB
f 2 FC
f 3FA   n1F 
 
f 3FB  n2 F 
f 3 FC   n3 F 
(38)
(39)
and finally,
s
FA
 r  q   q 1
FC
(40)
Extensive properties
F1 x1F1  F2 x1F2  FA x1FA  FB x1FB  FC x1FC
(41)
F1 x2F1  F2 x2F2  FA x2FA  FB x2FB  FC x2FC
(42)
F1 x3F1  F2 x3F2  FA x3FA  FB x3FB  FC x3FC
(43)
x1FA  x2FA  x3FA  1
(44)
x1FB  x2FB  x3FB  1
(45)
x1FC  x2FC  x3FC  1
(46)
For this case we have seven equations eq. (37) and eqs. (41) to (46) write in terms of r , s ,
q , molar composition ( xi F , xi j ), feed flow streams( Fk k  1, 2 ) , and product flow streams
k
( Fj
F
j  A, B, C ). Calculating equations – variables: 7-7, hence this equation system can be
solved. Expressing eqs. (41) to (46) as vector and matrix arrangements, and solving,
75
 x3FA  1  x1FA   x2FA 
 FB     FB   FB 
 x3   1   x1    x2 
 x3FC  1  x1FC   x2FC 
 
   
 FA   x1 A
 F    x FA
 B  2
 FC   x3FA
F
x1FB
x2 FB
x3FB
x1FC 

x2 FC 
x3FC 
1
  x1F1
  F1
  x2
  x3F1


x1F2 
F




x2 F2   1  
F
x3F2   2  
(47)
(48)
and s is calculated with eq. (40). Comparing eq. (39) with eq. (48), there is a clear advantage to
solve the material equations if this is made with intensive properties, the matrix expression is
very simple unlike equation (48) that has inverse matrix calculation. We need to find the
composition and flows in each section of the complex column as shown in Figure 23 for future
requirements at the composition profiles calculation. In the first section S 1 (see Figure 24), we
can obtain the V1 and L1 by the following relationship,
V1  (r  1) FA
(49)
L1  rFA
(50)
Figure 24. Mass Balance at the top section S1 in Figure 23 (b).
76
We know that the composition of vapor and liquid streams are equal to the product FA . So, we
can get yiV1  xiL1  xi FA . In the fourth section, S4 (see Figure 25), the flow rates V4 and L4 can be
described as the following equations,
V4  sFC
(51)
L4   s  1 FC
(52)
Figure 25. Mass balance at the bottom section S4 in Figure 23 (b).
The molar fraction composition of V4 ( yi FC ) is at equilibrium with FC , then yi FC is calculated by,
yiV4  Ki xiFC
(53)
where K i is the phase equilibrium constant of component i, K i is a function of the VLE model.
For constant volatility,
Ki 
i
c

j 1
for ideal solution,
j
xj
(54)
77
Pi sat
Ki 
P
(55)
and for nonideal case,
Ki 
 i Pi sat
P
(56)
After we obtain the vapor composition, we will be able to compute the liquid flow composition
by the following mass balance equation.
L4 xiL4  V4 yiV4  F C xiFC
(57)
Under the assumption of the constant molar overflow, the internal flow rates of the four column
sections can be related by
L2  L1  qF1
(58)
V2  V1  1  q  F1
(59)
L3  L4  qF2
(60)
V3  V4  1  q  F2
(61)
xiL2  xi L3  xiFB
(62)
V3  V2 , yiV 3  yiV 2
(63)
L2  L3  FB
(64)
The vapor composition of V2 is given by,
 yVA2  1 
 x AL1 
 x AL2 
 x AF1  

V

L

L

F



 V2 
1  L1 
2  L2 
1  F1  
 1

 yB  V2 
 xB 
 xB 
 xB  
(65)
78
2
yCV2  1   yiV2
(66)
i 1
Table III collects all specifications to find the vapor and liquid flows as well as vapor and liquid
compositions of the complex column shown in Figure 23.
Table III. Definitions of vapor flow, liquid flow, vapor and liquid compositions for all four
sections of a complex column
Section
VT
LT
yi ,T
xi ,T
1
V1  (r  1) FA
L1  rFA
xi FA
xi FA
First tray
composition
yi ,1  xi FA
2
V2  V1  (1  q) F1
L2  L1  qF1
yi ,V2
xi , L2
xi ,1  xi FB
3
V3  V4  (1  q) F2
L3  L4  qF2
yi ,V3
xi FB
yi ,1  xi FB
4
V4  sFC
L4  ( s  1) FC
yi ,V4
xi , L4
xi ,1  xi FC
2.2.2 Continuous Model for Complex Columns
The continuous model adapted to complex column is called difference point equation
model (Huss and Westerberg 1996; Huss and Westerberg 1996), where the column section is
redefined as a length of column between points of addition or removal of mass and/or energy.
The continuous model has the same structure such as a simple column model, but only things
necessary to be modified are boundary specifications to be applicable in all complex column
sections. Assuming constant molar flow, we have
Vyn1  LxT  Lxn  VyT
(67)
79
The liquid composition, xn , expands as a Taylor series related to the liquid composition, xn 1 by,
xn  xn 1 
dxn 1
n , where n   n  1  n  1. Rearranging equation (67) and defining the flow
dn
rate difference point   V  L (Holland et al. 2004), the reflux ratio difference point
dxi 
1 
1
 1 
 X i  xi 
  xi  yi  
dn  R 
R
(68)
The definition of  , R , and X i for any stage j is,
  Vj  Lj
2.3
R 
Lj

V j y i j  L j xi j
V
X i 
L

(69)
Temperature Collocated Continuous Profile Equation
Unfortunately, for higher space of multicomponent mixtures, a better model is needed for
a fully computational computer synthesis. We have discovered the thermodynamic
transformation of the column stage number, n, expressed in eq. (68) into temperature, T, as
shown in eq. (70). This temperature is equal to the tray temperature at each stage but also
corresponds to its bubble point, since vapor and liquid streams are in phase equilibrium. The
thermodynamically motivated transformation offers massive size reduction in the column profile
computations, eliminates singularities in the convectional composition profiles near stationary
and saddle pinch points, extends to any dimension of multicomponent mixture with any phase
equilibrium relationship as constant volatility, ideal, and nonideal mixtures, applies to both sharp
and sloppy splits; and is a key feature in a rigorous feasibility criterion based on minimum
bubble point distance functions which simplifies the search. A detailed derivation can be found
80
elsewhere (Zhang and Linninger 2004), including the nonideal and azeotropic versions of
equations.
c


1 
1
   1 
X j  x j 

xj  yj  
T
R 
R


K j
 T
x j
j 1
xj
 
 
1 
1
xj  yj  
X j  x j  K j 


  1 


R 
R
j 1 
 
 
c
(70)
In contrast to the restrictions of the conventional design methodologies, our approach
rigorously computes all liquid and vapor profiles resulting from the predicted operating
conditions such as reflux and reboil ratios, total tray numbers and the location of product
removal streams for all column sections, and a stream table of all intermediate and final products.
Once compositions profiles of the entire network have been calculated as a function of
temperature, the column height can easily be recovered by back transformation from temperature
to stage number as in eq. (71). Integration of the column profiles from the top to the bottom stage
gives the number of equilibrium stages, n. The number of stages is needed in capital cost
estimation
c
dn n x j


dT x j T
2.4
K j
 T
j 1
xj
 
 
1 
1
xj  yj  
X j  x j  K j 


  1 


R 
R
j 1 
 
 
c
(71)
Computation of Composition Profiles
The procedure to compute the composition profiles of all columns in a network is solving
eq. (70) as described in Table IV for each column section displayed in Figure 23. For each
column, the first step consists in specifying the product purities for a simple or complex column
81
Table IV. Pseudo-code implementation to compute the liquid composition profiles using
orthogonal collocation on finite elements for each column.
1.
2.
Specify the product purities, and design parameters such as reflux or reboil ratio
Calculate all product flow rates, internal vapor and liquid flow, and stationary points (Fout,
xj, yj, L, V)
3. Compute initial integration temperature, T
4. Calculate the column profiles parameters, R∆, X∆j, ∆
5. Compute the stable pinch point, xpinch j, Tpinch in each column section
6. Check if there is attainable temperature window between any pair of equivalent rectifying
and stripping column sections
7. IF the shortcut criterion is infeasible, THEN the column specifications are infeasible. Stop
8. Specify the placements and total number of elements and nodes
9. Create the temperature vector where each component of this vector denotes the
temperature value
N
d
T Te
10. Find the elements of the derivative matrix A m T  
lm T  lm (T )   e n e
dT
n 1 Tm  Tn
n m
11. Find g with equations  A  x    g 
12. Evaluate F(x)   A
 x   g 
13. Repeat steps 10, 11, 12, and 13 while F(x)  
14. If F(x en ,c )   , output the column section profiles
in term of molar composition, xPj, or fraction recovery, fPj, and design parameters such as reflux
ratio. Step 2 calculates all product flow rates by global mass balances in each column of the
network; this step is a linear algebraic problem. Steps 3 involves the computation of the
stationary nodes in each column section by the implementation of consecutive mass balances in
each section and bubble point temperature calculation; these stationary nodes are defined by the
values of xj, yj, T, L, and V. The calculation of the column profile parameters, R∆, X∆j, and ∆
is given in step 4. The stable pinch point of each section, xpinch j, Tpinch, is computed in step 5 to be
used in the next two actions. Steps 6 and 7 inspect a shortcut feasibility test based in the
existence of attainable temperature window (Zhang and Linninger 2004) between any pair of
82
equivalent rectifying and stripping column sections. If specifications pass the shortcut feasibility
test, steps 8 to 14 are executed to compute all liquid composition profiles of all species in each
section by solving eq. (70) between the stationary node, xj, T, and stable pinch, xpinch j, Tpinch,
using orthogonal collocation on finite element in each column section.
For regular rectifying sections, the generalized reflux becomes the classical reflux ratio,
and the difference point is equal to the distillate composition. In equivalent rectifying sections,
the difference point equation does not correspond to a product leaving the section. Yet, the
component balances for an equivalent rectifying section closes the component balances with the
difference point in the same way the distillate determines the operating line equations. For the
computation of a complex column section, we choose an intermediate product, whose
composition constitutes a starting point of the composition profiles. Accordingly, the
composition profiles described by eq. (70) for any two complex column sections emerge from
the respective intermediate product node. Similar considerations apply to equivalent stripping
sections.
2.4.1 Computation of Nonideal Mixture Model
In eq.(70), the equilibrium between vapor and liquid streams existing in equilibrium
stages is described by the equilibrium constant, Kj. For ideal mixtures, the equilibrium constant
for each species is given by the ratio of the pure component vapor pressure over the total
pressure, P. Accordingly, the K value changes with temperature along the column profiles. The
vapor pressures as a function of temperature are typically modeled by the well-known Antoine
equation. Nonideal phase equilibrium expressions including activity coefficients with implicit
dependency of composition and temperature are described by eqs. (72) and (73).
83
The nonideal vapor-liquid equilibrium model uses liquid phase activity coefficients. The
vapor phase fugacities were neglected to avoid the introduction of too many empirical
parameters. The temperature dependent binary parameters can be obtained using a commercial
database such as the Aspen properties tool (AspenTech). Analytical derivatives of the activities
coefficients with respect to composition and temperature are computed by general expressions
for several nonideal models (Taylor and Krishna 1993). To find the liquid composition profile,
the temperature continuous model equation in eq. (70), is solved for all sections belonging to
each column of any complex network following the procedure described in Table IV.
Kj 
PjV  j
P
v
 c  dx   
1  dPj
 
 j  Pjv   j i  j  
dT
P  dT
 i 1 xi dT T  
dK j
(72)
(73)
Eq. (70) can be solved implicitly after substituting eq. (73) resulting in a system of
ordinary differential equations. Using the temperature collocation for entire network requires no
modification of the minimum BPD and global feasibility test, since the liquid profiles are again
polynomial with piecewise elemental support whether the solution is ideal or nonideal.
2.5
Robust Pinch Point Computation
It has been demonstrated that design techniques such as the rectification body method,
column profile map, and the temperature collocation method all require the exact location of all
pinch points. Even points that lay outside physically realizable space have to be considered as
they determine the path of liquid composition profiles and delineate areas of reachable
compositions for the final design. This section is discussed for the efficient and rigorous location
84
of all Pinch Points, for the three VLE models; constant relative volatility, ideal mixture and
nonideal mixture models. We demonstrate here that pinch points for this model which is the
simplest of the three phase equilibrium models may be solved with polynomial arithmetic.
2.5.1 Pinch point location for constant relative volatilities
Pinch compositions can be found at the pinched condition, dx / dn  0 in the eq. (68).
R,i xi   R,i  1 yi  X .i  0
(74)
Since the constant volatility model relates the vapor and liquid phases as shown in eq. (26), eq.
(74) becomes by substituting yi
R ,i xi   R ,i  1
 i xi
 X  ,i  0
c
 x
j
j 1
(75)
j
Eq. (75)may be rewritten in terms of xi such that
c
xi 
X  ,i   j x j
j 1
R
c
 ,i
 1  i  R ,i   j x j
(76)
j 1
c
We now define a new parameter,     j x j , which can be interpreted as the mixture
j 1
vapor pressure, and is thus inversely related to the bubble point temperature This parameter is
also related to the Underwood roots, although Underwood never used this method for pinch point
location. Replacing  in Eq. (76) results in an expression for the Pinch Points as
xi 
R
 ,i
X  ,i 
 1  i  R ,i 
Using the fact the compositions sum to unity:
(77)
85
nc
 R
i 1
 ,i
X  ,i
 1  i  R ,i
1
(78)
By multiplying the denominator throughout in Eq. equation reference goes here, the equation
may be transformed into a polynomial function in terms of  :
anc nc  anc1 nc1  anc2 nc2 ...  a2 2  a1  a0  0
(79)
The polynomial coefficients a0 to anc are completely determined by specifying the
generalized reflux ratio RΔ and XΔ. The order of this transformed polynomial is equal to the
number of components in the system. This scalar equation is a simple, one-dimensional
polynomial, regardless of the amount components that are being considered. Once the
polynomial roots have been determined, the pinch compositions, xi, can be calculated from Eq.
(77).
The fact that the difference point equation reduces to a polynomial for a constant relative
volatility system, makes it entirely possible for the solutions to have imaginary parts. Thus,
although for a three component system there are always three solutions, only one of these may be
real as the other two are conjugate complex pairs. Even though two pinch points are complex, the
column profile with XT and YT as the entering composition has no odd or unrealistic behavior.
However, complex solutions are only found in systems where the difference point lies outside
the composition triangle (Tapp et al. 2004). A ternary column profile map with only one real root
is shown in Figure 26 for constant relative volatilities (A=5, B=2, C=1).
86
1.5
Real Solution
1
xC
0.5
X∆=[1.2,-0.1,-0.1]
XT
YT
0
XT=[0.4,0.22,0.38]
XΔ
YT=[0.9,0.02,0.08]
-0.5
-0.5
0
0.5
1
xA
1.5
R∆=0.6, L=1, V=2.67
Figure 26. A column profile map for a mixture with constant relative volatilities of ( A=5, B=2,
C=1) with only one real solution at a reflux of 0.6 and XΔ= [1.2, -0.1, -0.1]. The thick black
trajectory indicates a realizable liquid composition profile of a typical column section.
2.5.2 Pinch Point Location for for Ideal Solutions
Ideal solutions that adhere to Raoult’s law already elevate the complexity for solving the
pinch point equations due to the presence of an exponential term. Although a single solution may
be found with a basic numerical technique such as the Newton method, this firstly requires a
good initial guess of both composition and temperature. Secondly, even a good initial guess does
not guarantee convergence to the specific root that one wishes to find. Again using the bubble
point temperature will allow reducing the search for all pinch points to a single non-linear
equation. The difference point equation, for pinched conditions, may be rewritten in terms of the
liquid composition, such that
xi 
X  ,i
( R  1) Ki (T )  R
(80)
where Ki is the equilibrium constant defined as Ki (T )  Pi ,sat / P , and T is the bubble point
temperature. Again exploiting the fact that the sum of compositions is unity, the pinched
87
difference point equation may then be rewritten in the following form for an nc-component
mixture:
nc
X  ,i
i 1
( R  1) Ki (T )  R
f (T )  
1
(81)
Even though this function does not reduce to a polynomial function like in the case of
constant volatility system, this single non-linear equation, f(T), is significantly simpler to solve
than the original pinched difference point equation system in Eq. (74) because it only depends on
temperature and not composition (since the equilibrium constant Ki is a function of temperatures
only). Once all zeros are obtained in terms of the Bubble Point Temperature, the Pinch Point
compositions can be found by using Eq. (80). To better elucidate this scheme, in Figure 27 plots
Eq. (81) for a particular case of a quaternary mixture of pentane/hexane/heptane/octane, for a
specific reflux of 5 and a difference point (XΔ) of [0.9, 0.05, 0.03, 0.02].
Figure 27. A one dimensional pinch point temperature search showing all pinch temperatures
where f(T)=0 for XΔ=[0.9, 0.05, 0.03, 0.02] and RΔ=5, accompanied by a table showing all pinch
temperatures and compositions. The pinch temperatures indicated on the graph (T1-T4) are the
bubble point temperatures of the four pinch points.
88
2.5.3 One-dimensional deflation method for Pinch Point Calculation
Although the one-dimensional equation transformation described in the preceding
discussion significantly reduces the numerical complexity of the problem, a numerical technique
such as the Newton Method may still not find all desired zeros, especially if two zeros are
located nearby one another. To this end, a novel deflation method is proposed. By starting with
an initial guess for temperature greater than the maximum pinched temperature, convergence to
the highest pinch temperature is guaranteed. In order to locate all subsequent zeros quickly and
efficiently without convergence problems, we modify the original pinched Difference Point
Equation, f(T) in Eq. (81), to create a surrogate function, g(T), preventing the algorithm from
converging to a previous zero. The surrogate function, g(T), features previously found zeros in
the denominator, forcing it to near infinity near a previously found zero and thus preventing
repeated convergence to this specific solution. The sequence of surrogate equations, gnr(T), are
given in
g nr (T ) 
f (T )
nr
 (T  T
j 1
P, j
)
(82)
where nr is the number of zeros previously found and TP,j is a previously found pinch
temperature. This equation may then be solved in a relatively straightforward manner with a
simple numerical technique such as Newton-Raphson. Moreover, this method guarantees
convergence to all zeros as Eq. (82)nears infinity near previously found solutions. An illustration
of the systematic elimination of previously found zeros is shown in Figure 28. Figure 28 (a)
shows the surrogate function, g1(T), versus bubble point temperature after the highest bubble
89
point temperature has been eliminated. All subsequent solutions may then also be obtained by
systematically dividing by the product of all previous zeros as shown in Figure 28 (b) and (c).
(a)
10
(b)
0.2
(c)
4
0
T3
T2
T1
T2
2
-0.2
g2(T)
g1(T)
T1
g3(T)x10-5
0
-10
-20
-0.8
-40
-1
T1
-4
-0.6
-30
0
-2
-0.4
-6
-8
30
40
50
60
70
T(oC)
80
90
100
-10
30
40
50
60
T(oC)
70
80
90 35.3
35.5
35.7
T(oC)
35.9
36.1
Figure 28. (a) Surrogate function g1(T) vs. bubble point temperature where the highest pinch
point has been eliminated, (b) surrogate function g2(T) vs. bubble point temperature where the
two highest pinch points have been eliminated, and (c) surrogate function g2(t) vs. bubble point
temperature where the three highest pinch points have been eliminated.
2.5.4 Pinch Point Location for for Nonideal Solutions
Preceding sections have shown that, even for simple phase equilibrium models, the
difference point equation requires some modification in order to efficiently find all solutions.
Due to the fact that temperature and composition are implicitly linked to one another by the
activity coefficients, γi(x,T) a one dimensional equation transformation is not possible. This
implicit relationship between temperature and composition may be seen in the NRTL activity
coefficient model given in Eqs. (83) to (85):
  x j ji G ji
x j Gij

 i ( x, T )  exp  j

xk Gki
j  xk Gkj
 
k
k
 ij  aij 

 xm mj Gmj
  m
 ij
k xk Gkj







(83)
bij
T (K )
(84)
90
Gij  exp(cij ij )
(85)
Other numerical techniques such as Interval and Continuation methods have been applied
in an attempt to solve these non-linear equations. The interval method was shown to be a quick
and reliable tool for finding all azeotropes in both reactive (Maier et al. 2000) and non-reactive
systems (Maier et al. 1998). However, in these earlier works the search for azeotropic points was
confined to the mass balance triangle and global pinch point location was not pursued. Although
a Newton interval method may have potential to solve the Pinch Point problem, our trials with an
interval toolbox in MATLAB were not successful (Rump 2009). This is mainly due to the fact
that Newton interval method requires the entire interval to be smooth and differentiable, and the
multidimensional, highly non-linear nature of the NRTL equation is particularly difficult for
current interval computation techniques. Furthermore, even when confining the search interval to
differentiable regions, this method may not succeed in finding solutions that are located near one
another. Continuation methods, on the other hand, have the drawback that it requires pre-existing
knowledge of number and location of stationary points in the entire composition space and can
also be time consuming as an entire reflux locus has to be generated for a specified XΔ.
2.6
Piecewise polynomial approximation
This section describes the global search for the minimum bubble point distance function
for two adjacent column section profiles with polynomial arithmetic. The solution of eq. (70)by
finite element on orthogonal collocation generates an approximate analytical solution represented
by a set of polynomial temperature functions for each composition on each column section. Eqs.
91
(86) and (87) represent composition profiles for any pair of adjacent column sections r and s of a
multi-component mixture parameterized in terms of temperature T.
Xr T    g A T  , g B T  , gC T  , g D T 
Xs T    f A T  , f B T  , fC T  , f D T 
T
T
(86)
(87)
Then, a scalar function representing the Euclidean difference between any pair of
adjacent rectifying and stripping equivalent sections, BPD=||Xr(T)-Xs(T)|| evaluated at the same
temperature value defines the bubble point distance function, BPD. Hence, the BPD is also a
piecewise polynomial function in T of order equal to two times the highest degree of the two
polynomial orders. The BPD is a continuous scalar field in terms of temperature as in (88),
BPD T   X r T   X s T 
(88)
Using eq.(88), it is possible to compute continuous values of the BPD for any temperature
over which the pair of adjacent profiles exists simultaneously. This range is determined by the
temperature window given by the bubble point temperature of the first stage and the stationary
pinch point. Finally, the minimum bubble point distance, BPD, between two adjacent column
sections is the solution of a global optimization problem, (89)
  min BPD T 
T
(89)
Now, suppose N nodes are used per element to discretize the composition solution in the
temperature domain. In that case, the solution of the nonlinear profile equations are polynomial
functions gj(T) and fj(T) of order N-1 after implementation of eq. (88). Then, the BPD function is
a piecewise polynomial function of order 2(N-1) given as,
92
BPD(T )  a1T 2( N 1)  a2T 2( N 1)1 
 a2( N 1)1
(90)
To find its global minimum for a fixed design variable set and product requirements, the
optimality condition of eq. (90) generates another piecewise polynomial of one order less as
shown in eq. (91),
d  BPD T  
dT
  2( N  1) a1T 2( N 1)1  2( N  1)  1 a2T 2( N 1)2 
 a2( N 1)  0
(91)
The Newton Horner scheme (Zhang and Linninger 2006) with synthetic division can be
used to solve eq. (91) for all positive real zeros. The arrangement with the smallest BPD is the
desired minBPD. Note that polynomial arithmetic applies to solve the minimum bubble point
distance criterion globally even if the phase equilibrium models are non-linear.
2.7
Inverse Design Problem
In general, design of distillation columns can be performed by forward performance
simulation. It computes product distribution, vapor-liquid equilibrium relationship, liquid phase
enthalpies, and other outcomes relevant to column performance once we set thermodynamic
equilibrium model with given or priori calculated design specifications (e.g. feed, total column
tray, feed tray, reflux ratio, and etc). This approach requires trial-and-error adjustments of
operating conditions and column design specification, but it cannot guarantee desired product
purities because they are just results of the forward computation with input design specifications.
To obtain design specification for the specific product purity, it needs supplement procedure to
determine entire design and operating parameters or trial-and error methods. Our design
approach, on the other hand, starts from the designed product purities (both top and bottom
products) and inversely goes to find design parameters which was initially guessed from classical
93
performance design approaches (See the flowchart in Figure 29). Specifying product purities,
both column profiles of rectifying section (from top product) and stripping section (from bottom
product) is computed following methods described in sections 0 through 2.6.
For a feasible design these profiles should satisfy the following requirements. (i) the
search for configurations should incorporate all thermodynamically admissible combinations of
simple and complex column equipment to achieve any type of product distribution such as in
prefractionating, sharp, nonsharp, and sloppy splits; (ii) liquid composition profile equations
should admit all phase equilibrium relationships for ideal, nonideal, as well as azeotropic
mixtures; and (iii) column profiles must exactly intersect in feasible designs, while infeasible
design specifications are characterized by a gap in at least one pair of composition profiles in
adjacent column sections. Thermodynamically impossible specifications lead to nonzero bubble
point distances (Zhang and Linninger 2004). These requirements have been incorporated in an
inverse design methodology known as temperature collocation with rigorous computation of the
minimum bubble point distance (BPD).
2.7.1 Bubble Point Distance (BPD) for Feasible Design
After computing the liquid composition profiles of all column sections, the minimum
BPD is defined as the globally minimum distance between each pair of adjacent column sections.
In each section, the continuous temperature column profile equations are a set of polynomial
functions for each composition in terms of bubble point temperature as independent variable.
Since all composition profiles, xj(T), are approximated by polynomial functions, the minimum
BPD is found easily using polynomial arithmetic. A complex column k is feasible, if and only if
the sum of all minimum profile distances of any pair of equivalent rectifying, r, and stripping, s,
94
column sections is within a small -tolerance of zero, as in eq. (92), in which the sum of all its
column section bubble point distances is still a scalar. An entire separation network is feasible, if
and only if all its simple or complex columns, k, are feasible. The network feasibility is given in
eq. (93).
 (k )   min BPD(T )  1
(92)
K
 (k )    (k )   2
(93)
k 1
The procedure for meeting feasible design specifications for a single column is given in
Table IV. The information flow for the entire flowsheet is depicted in Figure 29. Product purities
of the final products are fixed in stage 1 of Figure 29. These purities can be specified in terms of
product compositions or fractional recoveries. The intermediate product purities and operating
conditions are at this point still to be determined so that the network becomes feasible and
optimal. In stage 2, initial guesses for these internal degrees of freedom such as intermediate
product compositions as well as reboil or reflux ratios of columns are specified. These choices
combined with suitable component balances allow the determination of all internal and exiting
process streams in the network. With all input and output streams of the network fixed, the BPD
criteria are easily computed for all pairs of column sections in the network. To find a feasible
network configuration, the minimum BPD algorithm finds the minimum distance between each
pair of adjacent column sections of the entire network in stage 3.
If all minimum distances are inside suitable tolerances close to zero, the entire network is
feasible. However, evaluating the BPD criterion for the initial guesses is not sufficient to
conclude a design. Different choices of the internal degrees of freedom may also be feasible with
95
even lower energy demand. Accordingly, the inner design optimization loop (stages 2 to 3) is a
global search. In case of infeasible specifications, these internal degrees of freedom are updated
by a stochastic methodology or by the design engineer until a feasible design is approached.
More details about the implementation of stochastic optimization for network synthesis can be
found elsewhere (Zhang and Linninger 2004). The inverse design methodology based in
temperature collocation is executed consecutively until all intersected column profiles are
resolved.
96
Stage 1
Set global feed composition, feed flowrate, and all
final product purities (xf, F, and f)
Stage 2
Set initial guess of intermediate products purities and
operating conditions for all columns in the network
(f, and r or s)
COmpute product flowrates, internal vapor and liquid
flowrates, and stationary points of all column sections
by species balances (Fout, V, L, xΔ, yΔ)
Specify different
operating condition
or intermediate
product purities
Calculate column profiles parameters and pinch points
(Δ, RΔ, and XΔ , xpinch, Tpinch)
Satisfy shortcut
feasibility criteria ?
Stage 3
no
yes
yes
Δ = V–L > 0 ?
for column sections
Equivalent rectifying section
integration from top to bottom
(TxΔ < Tpinch)
no
Equivalent stripping section
integration from bottom to top
(TxΔ > Tpinch)
BPD with polynomial arithmetic (eq. 2)
All BPD
satisfy tolerance?
no
yes
Feasible design
Figure 29. Detailed information flow diagram to find feasible design specifications for an entire
flowsheet. Stage 1 uses the feed and sets up the final product purities. In stage 2, initial guess of
intermediate products purities and operating conditions (reflux or reboil ratios) are specified.
Then product flowrate, internal vapor and liquid flowrate, stationary points and difference points
parameters of each column section are computed using species balances. Stage 3 performs BPD
criterion with obtained composition profiles by OCFE and polynomial arithmetic to assess
feasibility. The schematics on the right side depict the specifications for the computations of
each stage.
97
3
RESULTS OF DESIGN AND SYNTHESIS OF COMPLEX COLUMN
DISTILLATION NETWORKS
This section demonstrates the proposed complex column design methodologies using
inverse design procedure based on temperature collocation method for the separation of ternary
and quaternary mixtures in different ideal and nonideal mixture systems.
3.1
Ternary Mixture Separation with Complex Column Networks
For complex column system we used two arrangements, the first configuration has two
columns where the first is a single distillation column connected with a complex column similar
to Figure 30 to separate a constant volatility ternary mixture. All design and problem
specifications and vapor and liquid flows by section appear recollected in Table V. The column
profiles are calculated with Underwood equation model and Continuous model.
Figure 30. Complex column arrangement to separate a ternary mixture
98
Table V. Design and problem specifications and internal vapor and liquid flows.
Variables
Values

xF
[3.255; 1.809; 1]
[1/3; 1/3; 1/3]
r1
2.5
x F1
[0.6657; 0.3333; 0.0010]
x
F2
x FA
x FB
x FC
r2
V
L
[0.001; 0.3334; 0.6656]
[0.9800; 0.0199; 0.0001]
[0.0090; 0.9900; 0.0010]
[1e-9; 0.029999998; 0.970000001]
4.55
[1.8715; 1.8715; 1.8715; 1.8715]
[1.5343; 2.0342; 1.7147; 2.2148]
Figure 31 shows the column profiles using continuous model for a complex column network
shown in Figure 30. This is an example of feasible system because there is interception between
the profiles, where both have node pinches and the rectifying section has saddle pinch. The
saddle pinches correspond to a nearly infinite number of stages and a zone of nearly constant
composition in each column section (Doherty and Malone 2001). In contrast, at a pinch point the
compositions and temperatures are not function of the stage number. The two product streams
from the front column in Figure 30 are fed to the complex column which is the second column.
The complex system has three product stream and four profiles as in Figure 31 (b). Similar to the
simple column, the design and product specification produce a feasible complex distillation
column.
99
(a) Column Profiles for a Single Column
(b) Column Profiles for a Complex Column
Figure 31. Column profiles for a complex column network shown in Figure 30.
100
We also computed composition profiles with Underwood model for the same complex
column configuration of Figure 30. This figure show the characteristic profiles for a feasible
design and product specifications using the UW equation method, the profiles of the first column
are similar to the previously profiles calculated with continuous model, but this match is lose in
the complex column first profile. The overall column profile follows the path from the distillate
product, to the F1 feed tray, to the side-stream product ( FB ), to the F2 feed tray, and finally to
the bottom product. The extra bound in the graph, from the F1 feed tray to the F2 feed tray, is
due to the rectifying and stripping sections of the additional feed stream. The second
configuration consists of one complex column. Table VI shows the specification for this problem.
An important issue for mixture with three or more component is the difficulty to have problem
and design specifications that originate a feasible separation; this is an important motivation to
use optimization methods like genetic algorithm (Zhang and Linninger 2006) which is applied to
separation network systems.
Table VI. Design and problem specifications and internal vapor and liquid flows based on
Underwood method for the complex column configuration of Figure 30.
Variables
Values
Variables
Values

[3.255; 1.809; 1]
x FC
[0.00001; 0.0999; 0.9]
x F1
[0.6657; 0.3333; 0.001]
r
2.72
x F2
[0.001; 0.3333; 0.6657]
V
[1.4541; 1.4541; 1.4541; 1.4541]
x FA
[0.79; 0.209; 0.001]
L
[1.0632; 1.5632; 1.2907; 1.7907]
x FB
[0.09; 0.8; 0.11]
101
(a) Column Profiles for a Single Column
(b) Column Profiles for a Complex Column
Figure 32. Complex column arrangement to separate a ternary mixture
102
3.2
Quaternary Ideal Mixture Separation with Complex Column Networks
Finding design specifications which satisfy the aim of product purities as well as the
feasibility of the design is a necessity condition for the procedure of the realistic prefractionating
networks. A rigorous global energy minimization is another necessary task in the second stage of
inverse problem. Since the optimal design configuration is highly dependent on the mixture
condition, we used two different mixtures with two different feed mol-fractions. We had also a
case study for the nonideal mixture separation to show the robustness of our proposed method.
All mixture conditions used in case studies are listed as follows.

Alkane mixture of equimolar and different feed mol-fractions as (17.8%, 40.5%;
27.7% 14.0%) in Pentane, Hexane, Heptane, and Octane.

Aromatic and alkane mixture of equimolar and different feed mol-fractions as
(17.8%, 40.5%; 27.7% 14.0%) in Benzene, Toluene, Octane, and Nonane.

Nonideal mixture of equimolar fractions in Methanol, Ethanol, n-Propanol, and
Acetic acid (D).
As basic flowsheet configurations for quaternary mixture separation, we used 18
structures proposed by Agrawal (Agrawal 2003), in which includes 13 complex configurations
shown in Figure 33. When also admitting non-basic configurations without thermal coupling, 24
alternative configurations would be considered (Giridhar and Agrawal 2010; Giridhar and
Agrawal 2010).
103
D
Figure 33. 18 Separation flowsheets to separate a quaternary mixture. These configurations
contain simple columns or basic complex columns.
3.2.1 Case study I – Equimolar Alkane Mixture Separation
The purification of a quaternary alkane mixture of equimolar fractions in Pentane,
Hexane, Heptane, and Octane with vapor-liquid data given in Table VII was studied. From the
18 structurally distinct flowsheets in Figure 33, we chose 9 flowsheets, in which 4 simple
column networks and 5 complex column networks, to discuss detailed designs.
Table VII. Antoine parameters of the alkane mixture of Pentane (A), Hexane (B), Heptane (C),
and Octane (D) used for the calculation of the partial vapor pressure.
Components ( x A , xB , xC , xD )
[Pentane, Hexane, Heptane, Octane]
Antoine Parameter A
[6.85221; 6.8702; 6.8938; 6.9094]
Antoine Parameter B
[1064.630; 1168.7; 1264.4; 1349.8]
Antoine Parameter C
[232.000; 224.2100; 216.6360; 209.3850]
104
The selected flowsheets are depicted in Figure 34. There has a continuous feed stream of
100 kmol/h (27.78 mol/s) into products with purity of at least 99%. For each network, the
minimum BPD algorithm was used to determine feasible network specifications that met the
final purity constraints. Feasible column profiles were determined using temperature collocation
method as in section 0 and design specification were adjusted until feasibility and maximum
energy efficiency were achieved with the network BPD criterion in inequality eqs. (92) and (93).
NW1
NW4
NW7
NW2
NW5
NW8
NW3
NW6
NW9
Figure 34. 9 optimal basic network configurations to separate a quaternary mixture. These
configurations have an optimal tradeoff between capital and operating costs, and minimum
energy requirement.
105
The separation network configuration, NW1, has three columns. The first column operated
as a prefractionator renders sloppy splits and no final products. The second column is a ternary
simple column, which yields an indirect split with 99% purity of octane, component D, in the
bottom product. The third column consists of four complex column sections generating three
product streams at once, top-middle-bottom products. Figure 35 shows the composition profiles
of the entire network with precise intersection between all pairs of adjacent column sections in
each distillation tower.
P1
(a) flowsheet
F3
B
(b) Column C1
P2
F1
F4
P3
F2
F2
P4
F1
F3
D
C
A
(c) Column C2
B
(d) Column C3
B
P2
F2
F4
F4
D
D
F3
P4
C
C
A
P3
P1
A
Figure 35. Frame (a) shows the flowsheet with feeds (F1-F4) and products (P1-P4) for the
separation of a quaternary mixture. Frames (b) to (d) represent the tetrahedral composition
profiles of columns C1–C3 for the four species A, B, C, and D. The blue trajectories stand for
rectifying sections and red trajectories for stripping section.
For each column section, purities and refluxes had to be suitably adjusted to satisfy the network
BPD feasibility criterion. Global BPD minimization was performed with a stochastic global
106
search (Zhang and Linninger 2006) typically using the third component molar fraction in the
distillate product, x3,D and the reflux ratio, r, as the remaining search variable. Purities of internal
product cuts such as the distillate and bottoms for the prefractionator column were chosen to
lower the total energy consumption.
Note that traditional design methods assuming sharp splits are not suitable for the
synthesis of the heat integrated complex networks of Figure 35. Also the prefractionator is not
pinched so that minimum reflux or pinch alignment methods are inadequate for its design.
All nine network configurations shown in Figure 34 use quaternary columns as direct,
indirect, or prefractionator separation units were solved with the temperature collocation method.
The initial feed is a saturated liquid; the feeds to subsequent columns include saturated vapor for
lower energy consumption. Energy efficiency for each flowsheet was assessed by computing the
total vapor flow rate defined as the sum of the individual column vapor flows. Although total
vapor rate for the separation network does not account for different temperature levels required
in second law analyses, this parameter serves well as a first approximation of total expected
energy cost. The total vapor duty required for each of the nine configurations is given in Table
VIII. Table IX summarizes the design aspects and configurations of all nine networks computed
by our inverse design methodology. The synthesis rigorously synthesized networks leading to the
same final product purities.
The detailed network specifications in Table IX are used to initialize an equivalent
AspenPlus flowsheet. The Aspen performance calculation converged typically within a couple of
iterations in a few CPU seconds, thus confirming essential agreement between the two
computational results. The predicted column profiles of columns in networks NW 1, NW2 , and
107
Table VIII. Vapor rate for each column and total vapor rate of each designed network.
NW1
NW2
NW3
NW4
NW5
NW6
NW7
NW8
NW9
Col. I (kmol/h)
74.99
79.99
186.87
114.81
186.87
156.56
156.56
478.08
186.87
Col. II (kmol/h)
78.07
57.36
46.73
63.34
60.60
66.44
281.14
63.70
472.62
Col. III (kmol/h)
70.34
64.66
48.99
86.30
62.68
76.09
62.681
72.50
63.27
Total (kmol/h)
223.41
202.02
282.59
264.45
310.15
299.09
500.38
614.28
722.76
Table IX. Column design aspects obtained by solving the inverse design problem of nine column
configurations to separate a quaternary mixture using the minimum BPD design algorithm. In
this table, continuous length were converted to actual tray using eq. (71); remark were rounded
to next integer for better comparison with Aspen simulation.
Column I
Column II
Column III
NW1
NW2
NW3
NW4
NW5
NW6
NW7
NW8
NW9
Total stages
34
35
44
43
44
95
95
20
44
Feed stage
21
9
37
40
37
39
39
10
37
reflux
0.8
2.17
1.5
0.67
1.5
5.2
5.2
8.6
1.5
Total stages
70
31
39
40
44
43
28
27
22
Feed stage
39
9
16
12
15
18
20
11
16
reflux
1.36
1.09
0.498
1.51
1.4
1.2
4.68
1.548
8.5
Total stages
40
50
51
37
36
28
30
38
31
9
6
10
11
13
11
13
15
13
33
34
38
31
1.8
1.6
0.95
1.42
1.53
3.157
1.532
1.889
1.530
Feed stage
reflux
108
NW3 from our temperature collocation method are compared to the Aspen profiles in Figure 36,
Figure 37, and Figure 38 respectively. The graphs show the evolution of the liquid composition
profiles as a function of bubble temperature, Tstage, as well as in terms of stage number, NStage.
The mean error between the profiles obtained via the inverse and the rating problem is less than a
few percent. The difference is larger in the middle sections of the columns, but virtually the same
purities are obtained at the product nodes.
In the case studies reported here, the differences between the Aspen predictions and our
profile computations did not exceed 5%. The deviation at the products nodes was typically less
of 0.1%. Such small concentration differences are below the accuracy of standard concentration
measurement in the industrial practice, but may be significant in ultrapure gas separations. These
cases can be dealt with by a suitable tightening of the bubble point tolerances. We found the
agreement between our inverse designs with AspenPlus simulations remarkable, despite
qualitative differences between the continuous generalized profile equations used in temperature
collocation versus the tray-by-tray computations in flowsheet simulators. Six other network
configurations were also validated. We only show the composition profiles from temperature
collocation results of networks NW4-NW9 in Figure 34 for the sake of brevity. A critical
comparison confirmed the excellent agreement between our synthesis results and the precise
Aspen validation.
109
B
1
1
0.8
0.8
0.6
0.6
0.6
0.4
0.2
0
1
0.4
C
A
Xi
Xi
Xi
0.6
1
0.8
0
5
10
15
20
Nstage
25
C1
0.2
D
Xi
1
0.8
30
50
60
70
80
o
Tstage, C
90
0.4
0.4
0.2
0.2
0
1
100
0
5
10
15
20
Nstage
25
3032
1
1
1
0.8
0.8
0.8
0.8
0.6
0.6
C2
Xi
Xi
0.4
0.4
D
Xi
0.6
C
0.6
0.2
0.2
50
60
70
80
o
Tstage, C
90
100
Xi
1
0.4
0.4
0.2
0.2
B
0
20
30
40
Nstage
50
60 65
A
0.8
100
110
o
Tstage, C
Xi
B
0.2
0
1
C
10
15
20
Nstage
(a)
30
30
40
Nstage
50
60
0.8
0.8
0.8
0.6
0.6
0.6
C3
0.2
25
20
1
0
5
0
10
1
0.4
0.4
0
1
120 125
Xi
0.6
90
1
40
50
60
70
o
Tstage, C
80
90 97
90
100
110
o
Tstage, C
120 125
Xi
10
1
Xi
0
1
0.4
0.4
0.2
0.2
0
1
0
5
10
(b)
15
20
Nstage
25
30
40
50
60
70
80
o
Tstage, C
(c)
90 97
(d)
A
AB
C O MPLEX1
SIM PLE1
B
ABC D
BC
SIM PLE2
BC D
C
D
(e) Schematic flowsheet of NW1
(f) Schematic flowsheet of NW1 from AspenPlus
Figure 36. Initialization of AspenPlus by temperature collocation approach to validate network
NW1. Frames (a) and (b) show the composition profiles of all columns in the stage number and
temperature space respectively found by temperature collocation. Frames (c) and (d) describe the
profiles of the entire network in the stage number and temperature space respectively found by
AspenPlus. (remarks: component A: blue,--x--; component B: green,--●--; component C: pink,-■--; component D: red,--+--)
110
A
1
0.8
0.6
0.6
0.6
Xi
B
0.4
1
0.8
Xi
Xi
0.6
1
0.8
0.4
C1
C
0.2
0.2
Xi
1
0.8
0.4
0.4
0.2
0.2
D
0
25
30
0.8
C
50
60
70
o
Tstage, C
80
0.2
0
1
10
15
Nstage
B
0.8
0.6
0.6
B
25
C
C2
80
90
o
Tstage, C
100
0.2
0
5
10
15
20
Nstage
(a)
30
0.2
0.2
10
15
Nstage
20
25
29
1
0.8
0.8
0.8
0.6
0.6
0.6
C3
0
70
80
90
100
o
Tstage, C
(b)
(e) Schematic flowsheet of NW2
110
120
40
0
5
1
0.2
25
0.4
0
1
110
Xi
Xi
D
0.4
1
0.4
0.4
30 33
0.8
0.2
0.6
25
0.8
0.4
20
15
20
Nstage
0.8
0
5
10
1
D
1
0
5
1
0.6
Xi
0.4
0
1
90
Xi
0.6
40
1
50
60
70
o
Tstage, C
80
90
Xi
15
20
Nstage
80
90
o
Tstage, C
100
109
90
100
o
Tstage, C
110
123
Xi
10
Xi
5
1
Xi
0
1
0.4
0.4
0.2
0.2
0
1
5
10
15
20
Nstage
(c)
25
3032
0
70
80
(d)
(f) Schematic flowsheet of NW2 from AspenPlus
Figure 37. Initialization of AspenPlus by temperature collocation approach to validate network
NW2. Frames (a) and (b) show the composition profiles of all columns in the stage number and
temperature space respectively found by temperature collocation. Frames (c) and (d) describe the
profiles of the entire network in the stage number and temperature space respectively found by
AspenPlus. (remarks: component A: blue,--x--; component B: green,--●--; component C: pink,-■--; component D: red,--+--)
111
1
1
1
0.8
0.8
0.8
0.6
0.6
0.6
D
0.4
0.4
C
0.2
B
0
10
C1
0.2
A
0
Xi
Xi
Xi
0.6
20
Nstage
30
Xi
1
0.8
40
60
80
100
o
Tstage, C
0.4
0.4
0.2
0.2
0
1
120
0
5
10
15
20 25
Nstage
30
35
4042
1
1
1
0.8
0.8
0.8
0.8
0.6
0.6
Xi
C2
0.4
C
0.2
0.2
0
0
10
20
Nstage
B
A
0.8
50
30
1
C
Xi
10
15 20 25
Nstage
30
35
0.8
0.6
0.6
0
(a)
0
5
0.8
0.2
30
0
1
80
0.6
0
25
0.2
0.8
0.2
15
20
Nstage
0.2
1
0.4
10
0.4
1
0.4
5
0.4
1
C3
Xi
0.6
60
70
o
Tstage, C
40
50
60
70
o
Tstage, C
(b)
(e) Schematic flowsheet of NW3
80
90
70
80
90 100
o
Tstage, C
110
120125
50
60
70
o
Tstage, C
50
60
70
o
Tstage, C
8082
Xi
A
Xi
0.6
B
Xi
0.4
Xi
0.6
60
Xi
1
0.4
0.4
0.2
0.2
0
1
0
5
10
15
20
Nstage
(c)
25
30
34
40
80
90 95
(d)
(f) Schematic flowsheet of NW3 from AspenPlus
Figure 38. Initialization of AspenPlus by temperature collocation approach to validate network
NW3. Frames (a) and (b) show the composition profiles of all columns in the stage number and
temperature space respectively found by temperature collocation. Frames (c) and (d) describe the
profiles of the entire network in the stage number and temperature space respectively found by
AspenPlus. (remarks: component A: blue,--x--; component B: green,--●--; component C: pink,-■--; component D: red,--+--)
`
112
Tstage, oC
Tstage, oC
95 105
Tstage, oC
56 78
85
0.8
92 107118
1
58 87
0.8
C
0.6
5
0.8
15
41 50
20 25
Nstage
30
Tstage, oC
55 65
7579
1
5
36
10
15 20 25 30
o
Nstage
Tstage,
C
Tstage, oC
41 46 56
61 71
35
7681
0.4
B
C
89
91
25
30
31
104
110
10
83
15 20
Nstage
Tstage,
98
25
30
oC
100
125
1
0.8
Xi
C
0.2
0
1
36
D
0.4
B
0.2
C
C2
Xi
Xi
C2
69
5
10
15
20 25 30 35
Nstage Tstage, oC
Tstage, oC
71 74 78 87 94
0
1
40
98
1
0.8
B
5
10
99
15
20 25 30 35 o40 43
Tstage, C
Nstage
Tstage, oC
100
104
110 119
125
0.8
B
C
0.6
Xi
C3
0.4
0.6
Xi
0.6
C3
0.4
C
0.4
C
0.2
10
15
20
Nstage
D
0.2
D
5
77
15
20
Nstage
o
Tstage, C
Tstage, oC
0.6
C
69
69
10
B
0.4
5
5
0.8
A
0.6
0
1
D
0
1
40
41
1
0.8
0.6
0
C
0.2
B
0
1
35 40 44
A
0.2
A
0.4
A
10
36
B
D
0.2
90
Xi
C1
0.4
A
0
1
Tstage, oC
60
72
78
58
0.6
D
B
36
0.8
C
Xi
Xi
C1
0.2
1
1
0.6
0.4
1
109 125
Xi
1
Tstage, oC
25
30
0
1
0.2
5
10 11
0
1
Nstage
5
10
14
Nstage
(a) Composition profiles from temperature collocation of NW4, NW5, and NW6
(b) Schematic flowsheet of NW4, NW5, and NW6
Figure 39. Composition profiles versus number of stages and temperature stage domains found
by temperature collocation for networks NW4, NW5, and NW6. Frame (a) shows the composition
profiles of all columns in the stage number and temperature space for NW 4 , NW5, and NW6.
Frame (b) contains the schematic flowsheet diagrams for all three networks. (remarks:
component A: blue,--x--; component B: green,--●--; component C: pink,--■--; component D:
red,--+--)
113
o
Tstage, C
36
48
8590
1
0.8
48
Tstage, oC
69 75
62
0.8
A
B
1
C
B
Xi
C1
0.4
C
C1
0.4
0.4
A
0.2
0.2
B
D
1
10
15
20
Nstage
A
0
1
25
31
o
Tstage, C
Tstage, oC
99 104 113 122 125
81 9194
5
10
Nstage
Tstage,
1
0.8
D
D
0.2
5
36
40
43
15
0
1
18
5
10
15
o
oC
47
Tstage, C
20 25
Nstage
Tstage,
50
65 68
1
oC
48 59 64
C
0.8
C
Xi
C2
0.4
D
0.4
B
B
A
0.2
5
10
15
Nstage
69
71
Tstage, oC
74
80 81 93
0.2
0
1
20
24
Tstage, oC
98
1
5
99
100
10
Nstage
Tstage, oC
104106
0
1
15
Tstage, oC
109119
125
1
0.8
B
Xi
C3
0.4
0.4
B
0.2
10
Nstage
15
17
0
1
51 65 68
A
D
C
15
19
Tstage, oC
0.6
Xi
Xi
C3
0.4
5
36
10
Nstage
Tstage, oC
40
47
0.8
0.6
0.2
5
C
0.6
0
1
98
0.6
Xi
Xi
C2
0.8
41
o
0.8
0.6
0.4
1
35
Tstage, C
70 75 85 93
A
0
1
30
B
0.6
0.2
109 125
0.6
C
0
1
Tstage, oC
95 105
58 87
0.8
0.6
Xi
0.6
87 97 102 109
Xi
1
Tstage, oC
60
72 78
0.2
5
10
Nstage
15
18
0
1
5
10
Nstage
15
17
(a) Composition profiles from temperature collocation of NW7, NW8, and NW9
(b) Schematic flowsheet of NW7, NW8, and NW9
Figure 40. Composition profiles versus number of stages and temperature stage domains found
by temperature collocation for networks NW7, NW8, and NW9. Frame (a) shows the composition
profiles of all columns in the stage number and temperature space for NW 7 , NW8, and NW9.
Frame (b) shows the schematic flowsheet for each network. (remarks: component A: blue,--x--;
component B: green,--●--; component C: pink,--■--; component D: red,--+--)
114
3.2.2 Case study II – Non-equimolar Alkane Mixture Separation
Unlike previous case study in section 3.2.1, we studied complex column separation of
alkane mixture of Pentane, Hexane, Heptane, and Octane with non-equimolar feed condition. By
changing feed condition with of [0.1780; 0.4046; 0.2771; 0.1403], we found the new design
parameters to separate this mixture as given in Table X.
Overall, each column for each network had a small variation of total number of stages but
different reflux ratio. These changes affect to the total energy consumption throughout the
network. For an example, the reflux ratio of the second column was increased from r=2.0 to 3.5
and the energy consumption for the last column could be reduced from r=2.46 down to 1.77 in
the case of NW1. On the other hand, the total numbers of stages in the second and third columns
had less variation down by 5 and 1 respectively. It is interesting to note that even though the
thermodynamics are completely different, the design methodology can be easily used to obtain in
detail all design specifications such as product streams, feed and product stages as before. In
addition, we can infer that design parameters must be determined highly depending on the
condition of feed into the distillation process, otherwise the efficiency of the process will be
attenuated.
A comparison of the energy consumption to separate this mixture between the two feed
scenarios shows that the vapor duties for both feed conditions were similar requiring about
200kmol/h in NW1. In cases of other networks, vapor duties are relatively different from
networks and feed conditions, but NW1 had a least vapor duty for the separation of this mixture
for both feed conditions. This result was shown in Figure 42.
115
Table X. Detailed network design specifications and operating conditions obtained by solving
the inverse design problem to separate a quaternary mixture of Pentane, Hexane, Heptane, and
Octane using two different molar compositions of equal molar FEED1 of [0.25; 0.25; 0.25; 0.25]
and FEED2 of [0.1780; 0.4046; 0.2771; 0.1403] in three different complex networks NW1, NW2,
and NW3 as in Figure 34.
NW1
FEED1
FEED2
FEED1
FEED2
reflux
21
1
43
43
100
37.78
62.22
0.500
8
36
1
21
40
40
37.78
62.22
25.00
25.05
49.95
2.000
22
1
40
40
100
38.21
61.79
0.620
6
27
1
20
35
35
38.21
61.79
17.94
40.30
41.76
3.5
18
1
21
21
100
52.21
47.79
1.800
12
28
1
15
30
30
52.21
47.79
49.85
24.95
25.20
3.900
21
1
23
23
100
78.84
21.16
0.600
13
33
1
16
35
35
78.84
21.16
57.94
27.98
14.08
2.700
Feed stage location
16
16
7
24
Product stage location
1
23
1
22
1
24
1
30
Total stages
23
22
24
30
Feed Flows (kmol/h)
49.95
41.76
49.85
57.94
Product Flows (kmol/h)
24.82
25.13
27.73
14.03
25.05
24.80
17.60
40.34
Reflux
2.460
1.766
1.457
2.543
17
1
21
21
100
50
50
0.700
6
22
1
16
23
23
50
50
24.86
28.35
46.79
4.100
8
29
1
20
39
39
28.35
46.79
25.18
24.96
25.00
3.800
18
1
20
20
100
56.97
43.03
1.200
7
33
1
19
36
36
56.97
43.03
17.85
43.08
39.07
2.7
11
29
1
16
33
33
43.08
39.07
40.26
27.74
14.15
1.900
Product stage location
Total stages
Feed Flows (kmol/h)
Product Flows (kmol/h)
reflux
Feed stage location
Product stage location
C2
Total stages
Feed Flows (kmol/h)
Product Flows (kmol/h)
C3
NW3
FEED2
Feed stage location
C1
NW2
FEED1
116
Figure 41. Liquid composition profiles on each column of the networks NW1, NW2, and NW3.
These feasible designs were found by temperature collocation method. The quaternary mixture
MIX2 constituted by Pentane (A), Hexane (B), Heptane (C), and Octane (D) is separated in each
one of the three networks with feed, FEED2 of [0.1780; 0.4046; 0.2771; 0.1403]. The
composition profiles on the top are corresponding to first column, second row to second column,
and bottom profiles to third column for each network. In the composition triangle, the labels, Fi,
refer to the feed tray compositions, not the feed compositions.
117
Figure 42. Total vapor flow rate required for the separation of alkane quaternary mixture of
Pentane, Hexane, Heptane, and Octane using two different molar compositions of equimola feed
of [0.25; 0.25; 0.25; 0.25] and non-equimolar feed of [0.1780; 0.4046; 0.2771; 0.1403] in three
different complex networks NW1, NW2, and NW3.
3.2.3 Case study III – Aromatic and Alkane Mixture Separation
For the robustness of the proposed design methodology and finding the energy efficient
design specification of specific petroleum-based mixture, we have considered the separation of a
mixture containing aromatic as well as alkane compounds; Benzene (A), Toluene (B), Octane (C),
and Nonane (D). Two different separation problems corresponding to feed compositions as
equimolar feed and non-equimolar feed of [0.178; 0.4046; 0.2771; 0.1403] were studied again.
The networks from NW1 to NW3, shown in Figure 34 were tested for the purpose of energy
efficiency of complex column network. The vapor liquid equilibrium model shown in Table VII
was used, and the aim of distillation was set capable of delivering the products with the same
99% purity requirement.
118
In result of temperature collocation method as in Table XI, the main parameter
modifications occur in the total trays and refluxes for the first and last columns. For the case of
non-equimolar feed condition, the first column runs at a 43% higher reflux than equimolar feed
codition. Yet, the third column needs 9 less stages with the reflux diminished by 28%. The
second column can keep similar reflux and total stages. Overall, NW1 requires 11% less energy
to separate FEED1 than FEED2. The investment cost, expressed roughly by considering tray
numbers only, suggests that FEED2 requires less capital cost.
In cases of NW2, it deploys a non-sharp split prefractionator C1 followed by a complex
column C2. Here, the component C (octane) is split sloppily. In contrast to NW1, the complex
column C2 products C and D (nonane) splits, while the lights A (benzene) and B (toluene) are
purified in column C3. The purities of the final products produced by NW 2 satisfy the
specifications as desired. Network 3 has two complex columns including several non-sharp
splits. It engages column C1 for a sloppy split of two components, B and C. The complex column
C2 produces only benzene in high purity. The complex column C3 finally completes the task by
purifying simultaneously toluene, octane, and nonane.
Composition profiles of both equimolar and non-equimolar feed cases for all three
networks are depicted in Figure 43 and Figure 44 respectively. Both feed condition cases have
marked differences in the intermediate streams as well as the detailed design specification for
each column. For example, the prefractionator purity leaving C 1 at the top is much richer in A,
but the concentration of B is 50% less than for the first mixture MIX1. Finally, the purity of the
terminal products is almost same, consistent with the high purity requirements (>99%). The final
119
molar flow rate of product B, C is larger for FEED2 as is expected given the different amounts in
inlet molar stream of FEED2.
Table XI. Detailed network design specifications and operating conditions obtained by solving
the inverse design problem to separate a quaternary mixture of Benzene, Toluene, Octane, and
Nonane using two different molar compositions of equal molar FEED1 of [0.25; 0.25; 0.25; 0.25]
and FEED2 of [0.1780; 0.4046; 0.2771; 0.1403] in three different complex networks NW1, NW2,
and NW3 as in Figure 34.
NW1
FEED1
FEED2
FEED1
FEED2
reflux
25
1
42
42
100
38.88
61.12
1.675
8
45
1
22
58
58
37.88
62.22
25
25.05
49.95
4.100
26
1
43
43
100
52.54
47.46
2.392
11
41
1
20
55
55
52.54
47.46
17.82
40.25
41.93
4.200
17
1
28
28
100
61.25
38.75
1.600
24
55
1
42
64
64
61.25
37.75
49.78
25.40
24.82
3.700
17
1
29
29
100
70.67
29.33
1.550
26
48
1
37
66
66
70.67
29.33
57.89
28.15
13.96
3.800
Feed stage location
10
17
14
9
Product stage location
1
34
1
25
1
20
1
20
Total stages
34
25
20
20
Feed Flows (kmol/h)
49.95
41.93
49.78
57.89
Product Flows (kmol/h)
24.82
25.13
27.90
14.03
25.03
24.75
17.60
40.29
25
26
1.928
3.308
21
1
32
32
100
44
56
0.950
7
32
1
14
35
35
44
56
25.21
43.32
31.47
3.000
14
59
1
45
63
63
43.32
31.47
25.37
24.29
25.13
3.5
22
1
34
34
100
56.06
43.94
1.500
13
35
1
24
47
47
56.06
43.94
17.97
46.63
35.40
2.900
16
51
1
40
61
61
46.63
35.40
40.49
27.61
13.93
2.7
Product stage location
Total stages
Feed Flows (kmol/h)
Product Flows (kmol/h)
reflux
Feed stage location
Product stage location
C2
Total stages
Feed Flows (kmol/h)
Product Flows (kmol/h)
C3
NW3
FEED2
Feed stage location
C1
NW2
FEED1
reflux
120
o
1
130
1
1
0.6
0.6
0.4
0.4
0.4
0.2
0.2
0.2
10
20
Nstage
30
0
40
5
10
15
Nstage
o
85
20
0
25
111
127
133
1
92
117
148
151
1
0.6
0.6
0.4
0.4
0.4
0.2
0.2
0.2
20
30
Nstage
40
0
50
10
20
o
126
30
40
Nstage
50
150
1
80
110
1
0.6
0.4
0.4
0.4
0.2
0.2
0.2
0
0
30
P1
F2
5
10
Nstage
F1
C1
15
20
A
P1
AB
ABCD
F4
P2
C2
F3
P3
BCD
F4
CD
10
C3
111
20
Nstage
30
Tstage, C
125
116
AB
ABCD
F1
C1
C2
0
10
20
A
148 151
30
40
Nstage
F2
ABCD
P2
P3
A
B
C
F1
C1
C2
F3
P2
F4
B
BC
C3
BCD
P3
C
F5
CD
P4
(b)
60
P1
CD
D
50
ABC
C3
ABC
F3
P4
(a)
F2
B
C
133 141
Xi
0.6
Xi
0.6
Xi
0.8
25
30
o
96
0.8
15
20
Nstage
25
Tstage, C
116
90
0
60
0.8
10
80
Tstage, C
134
5
15
20
Nstage
o
Tstage, C
1
10
Xi
0.6
Xi
0.8
Xi
0.8
10
131
o
Tstage, C
126
0.8
0
5
o
Tstage, C
80
Tstage, C
113
91
Xi
0.6
Xi
0.8
1
C3
140
0.8
0
C2
o
Tstage, C
119
96
0.8
Xi
C1
o
Tstage, C
113
88
D
P4
D
(c)
Figure 43. Liquid composition profiles on each column of the networks NW1, NW2, and NW3.
These feasible designs were found by temperature collocation method. The quaternary mixture
of Benzene, Toluene, Octane, and Nonane is separated in each one of the three networks with
equimolar feed, FEED1, of [0.25; 0.25; 0.25; 0.25]. The composition profiles on the top are
corresponding to first column, second row to second column, and bottom profiles to third column
for each network. In the composition triangle, the labels, Fi, refer to the feed tray compositions,
not the feed compositions.
121
o
1
129
1
1
0.6
0.6
0.4
0.4
0.4
0.2
0.2
0.2
0
0
10
20
Nstage
30
40
5
10
15
Nstage
o
80
102
111
20
0
25
130
1
98
117
151
1
0.4
0.4
0.4
0.2
0.2
0.2
30
Nstage
40
0
50
10
20
o
126
30
40
Nstage
50
136
150
1
110
1
0.4
0.4
0.4
0.2
0.2
0.2
15
Nstage
20
25
P1
F2
0
5
10
Nstage
F1
C1
15
20
A
P1
AB
ABCD
F4
P2
C2
F3
P3
BCD
F4
CD
Tstage, C
112
127
101
10
20
30
Nstage
C3
111
131
40
Tstage, C
126
113
AB
ABCD
F1
C1
C2
0
10
20
A
141
151
30
40
Nstage
F2
ABCD
P2
P3
A
B
C
F1
C1
C2
F3
P2
F4
B
BC
C3
BCD
P3
C
F5
CD
P4
(b)
60
P1
CD
D
50
ABC
C3
ABC
F3
P4
(a)
F2
B
C
30
Xi
0.6
Xi
0.8
0.6
Xi
0.8
10
25
o
Tstage, C
96
80
0.6
5
80
0
60
0.8
0
15
20
Nstage
o
Tstage, C
1
10
Xi
0.6
Xi
0.8
0.6
Xi
0.8
20
129
o
Tstage, C
125
146
0.6
10
5
o
126
0.8
0
Tstage, C
115
99
Xi
0.6
Xi
0.8
Tstage, C
C3
136
0.8
1
C2
o
Tstage, C
117
102
0.8
Xi
C1
o
Tstage, C
114
97
D
P4
D
(c)
Figure 44. Liquid composition profiles on each column of the networks NW1, NW2, and NW3.
These feasible designs were found by temperature collocation method. The quaternary mixture
of Benzene, Toluene, Octane, and Nonane is separated in each one of the three networks with
non-equimolar feed, FEED2 of [0.1780; 0.4046; 0.2771; 0.1403]. The composition profiles on
the top are corresponding to first column, second row to second column, and bottom profiles to
third column for each network. In the composition triangle, the labels, Fi, refer to the feed tray
compositions, not the feed compositions.
122
An important feature is that for FEED1, NW3 is more energy efficient than NW1. The
energy savings -without the claim for rigorous optimality- amount to 10%. From the energy
point of view, configuration NW3 is best to separate FEED1. Interestingly for FEED2, NW1 and
NW3 have almost the same energy requirement. We also infer from Figure 45 that NW2 seems to
be suboptimal for either feed composition from an energy consumption point of view. The
detailed comparisons of the structural alternatives for the feed compositions, FEED1 and FEED2,
suggest that different structural network solutions have a strong impact on the overall energy
consumption required for addressing a separation task. Moreover, our methodology seems
suitable for performing these detailed comparisons, as the purities are ensured for realistic
columns without dropping rigor in the calculations or the need to employ short-cut methods or
the restrictions of sharp splits. The networks are really comparable as they all solve with the high
precision identical separation tasks.
Figure 45. Total vapor flow rate required for the separation of aromatic and alkane quaternary
mixture of Benzene, Toluene, Octane, and Nonane using two different molar compositions of
equimola feed of [0.25; 0.25; 0.25; 0.25] and non-equimolar feed of [0.1780; 0.4046; 0.2771;
0.1403] in three different complex networks NW1, NW2, and NW3.
123
DISCUSSION
It was demonstrated that the temperature collocation algorithm combined with minimum
BPD can adjust the design specifications so that the product purity targets are met, if these
targets are at all attainable. Using the fractional recoveries as degrees of freedom, we enforced
global column species balances. Columns species balances would thus close provided that their
respective profiles pass the rigorous feasibility test. Suitable fractional recoveries can be selected
by the designer. The search space is usually small because of the specified high product purities.
Recoveries would typically range in the region of 0.99 and higher. Alternatively, a black box
search can be performed automatically using suitable global search algorithms such as the hybrid
algorithm described previously. Unfortunately, there are no obvious criteria to identify feasible
network recoveries a-priori, accordingly our semi-automatic design feasibility test needs to be
deployed repeatedly or a global search can be run as described elsewhere.
In the results of total vapor flow rates of each network, all products meet the high purity
separation targets as 99%. Network NW1 has among all structures the least total vapor rate for
the separation of MIX2 . With almost the same amount of energy - with less than 1% difference both feed compositions of FEED1 and FEED2 for MIX2 can be purified. Our synthesis method
precisely identified the necessary adjustments in terms of number of trays as well as locations of
feed and side-products. These detailed results could also permit a detailed capital cost assessment
combined with more explicit operating cost analysis. These types of detailed cost optimization
are left to the industrial practitioners, but are ideal follow-up steps intended to be performed after
the use of our synthesis methodology.
124
For mixture MIX1, NW3 has the lowest energy consumption. This configuration, NW3, can
again split both feed compositions with about 20% of the energy demand. NW1 and NW3 have
almost the same energy requirement for splitting FEED2 of mixture MIX1 with only 2.6% of
difference in terms of total vapor flow. NW2 has the worst energy performance in all cases,
because the two feeds specifications do not benefit from the pseudo-indirect separation occurring
in this particular network. It is also the only system, in which the equal molar mixture, MIX2FEED1, consumes more heat than the second feed composition, MIX2-FEED2.
Another interesting result is given by the total vapor consumption in NW3 to solve the
four different problem specifications. This network needed an average vapor rate of 313.5kmol/h
for all four cases with only 15% deviation from the average vapor consumption in the system
MIX1-FEED2. In this example, its performance may be explained by its intensive degree of
heat-integration, leading to almost constant heat duty for any separation task in our examples.
However, we recommend not to draw general conclusions from specific examples like the ones
we show here, but to consider the rigorous synthesis and analysis results for each particular case
as a function of the different thermodynamic systems with their specific feed compositions and
purity requirements.
The perfect test for the validity of the synthesis results would be to build the pilot or
process plants in the real world. Since this approach would lead to unreasonable cost, we
resorted to the next best validation based on commercial flow simulators which solve the
rigorous mass, equilibrium, summation, and heat equations (MESH). Even though the
temperature collocation approach approximates the discrete MESH equations by continuous
generalized profile equation, a remarkable agreement between our synthesis and the flowsheet
125
simulation occurred in all cases we investigated. The achieved purity targets in both methods
were practically identical, and the flowsheet simulation converged in a few iterations, in essence
to the same profiles as computed by our approach. Small differences in the middle of column
section trajectories obtained by temperature collocation and Aspen are due to model differences
like discrete versus continuous profiles as well as possible heat effects not accounted for in our
method, while Aspen does.
The inverse design methodology based in temperature collocation can also be used to
effectively initialize flowsheets for rigorous flowsheet simulation like AspenPlus. The agreement
between the converged AspenPlus product specifications and the specified product requirements
indicates that temperature collocation synthesis algorithm is practically equivalent to AspenPlus
results.
126
4
DISCUSSION
In this study, we aimed at novel methods for discovering energy-efficient distillation
configurations with computer-based approach. Distillative separation popular in the petroleum
industry deals with mixtures of tons of components, but systematic approach on quaternary
mixture still has not been fully accomplished and therefore it is necessary to find optimal
operating conditions and configuration to separate quaternary mixture or sometimes more
components mixture with the index of smaller energy consumption. Furthermore, the automatic
synthesis of multicomponent mixture separation especially for nonideal solution models remains
an unsolved problem for conventional design methodologies.
Thus far, conventional methods to design multicomponent mixtures often use
pseudocomponents by lumping components in the feed to reduce the size of the system. This
approach reduces the mixture to be separated into a set of pseudocomponents upon which
algorithms of minimum reflux conditions can be applied for sharp separations or related species
mixtures. However, lumping rules are limited to grouping components with similar volatilities.
Usually, nonkey components are lumped, leading to sharp split assumptions for nonkey
compounds which do not achieve realistic results. It also assumes that components’ molar
fractions in the lumped feed and product cuts are equal. Lumping criteria should not be applied
in cases such as the following: i. individual concentration in a pseudocomponent such as a highly
toxic material (e.g., sulfur compounds in crude oil like H2S43) should be purified, ii. the relative
volatility of the individual components to be lumped are very different, and iii. the amount of
127
species to be grouped in the pseudocomponent in the feed are not exactly the same as the internal
composition in the pseudocomponent in the product.
In these situations, rigorous and complete simulation of multicomponent mixtures was
advisible. Temperature collocation drastically reduced the problem size for separation design, so
that species lumping was not necessary. Moreover, the inverse design methodology based on
temperature collocation elucidated the topology of the thermodynamic search space in separation
problems to identify feasible operating conditions for multicomponent separation problems
including nonideal mixtures, prefractionating column configurations, and high-order
multicomponent mixtures.
Providing a number of examples of heat integrated complex networks using
prefractionators, we arrived at the important conclusion that the synthesis of optimal networks
fortunately cannot be performed with shortcut methods, performed with rule-of-thumb or
combinatorial separation heuristics, or addressed by intuition. For a complex separation network
solution to be optimal, the separation problem needs to be solved rigorously for the specific feed
compositions and the appropriate phase equilibrium model.
In addition, all design specifications obtained by temperature collocation were validated
by the commercial flowsheet simulation. Initializing AspenPlus flowsheet simulations converged
in a few iterations, achieving essentially the same product requirements without the need to
adjust the design specifications obtained by our inverse design methodology. The results show
that designs obtained by temperature collocation are equivalent to rigorous MESH solutions.
Alternatively, the design method can be considered a robust and rigorous process to initialize
AspenPlus flowsheet simulations.
128
Current study, however, may not been used for azeotropic separations. The synthesis of
azeotropic separation networks, however, requires a more complex logic, which this article does
aim to address. We have so far not attempted to solve for optimal designs, while enforcing
feasibility. This task would involve a global multilevel optimization with a high-level objective
function subject to structural (which flow sheet, which configuration) and parametric (operating
conditions, actual product purities) decisions. The feasibility criterion would have to be enforced
as a side condition; this sub-problem is already a global optimization problem.
129
5
CONCLUSION
The applicability and stringency of the minimum bubble point distance algorithm
combined with temperature collocation for the successful design of heat-integrated complex
separation networks with various feeds and volatility differences is studied. The minimum BPD
criterion was robust and efficient to characterize feasible network specifications guaranteeing
continuous liquid composition profiles connecting all products and feed stages in each column
section of complex flow sheets. Specifications leading to a nonzero bubble point distance are
guaranteed to be impossible to realize. The proposed feasibility test is therefore rigorous, subject
of course to the limits of floating point operations on a digital computer.
To separate four component petroleum mixtures, we simulated 18 different heatintegrated networks including simple column networks. Each design and synthesis problem was
rigorously solved to compute key design parameters. Our results suggest that it is not possible to
generalize that a particular network structure could be the best in terms of energy requirements
without studying in detail the separation task as a function of feed, the thermodynamic system,
and the feed composition of the mixture.
Finally, a remarkable agreement between the inverse design solutions and rigorous
forward flow sheet simulations was demonstrated. And possibly feasible designs were validated
by initialization of AspenPlus simulations which converged in a few iterations in excellent
agreement with the solution from the temperature collocation method. The converged product
specifications in the Aspen simulations match closely the final product purities computed by
130
inverse design. These results show that temperature collocation can be used to automate the
synthesis of complex separation networks.
This study suggests that the synthesis methodology can be adapted to address modern
process design problems such as computer-aided synthesis of energy-efficient separations,
detailed design of large-scale biorefineries with novel feedstocks, or new sustainable process
designs for reduction of greenhouse gas emissions.
131
CITED LITERATURE
Agrawal, R. (2003). "Synthesis of multicomponent distillation column configurations." AIChE
Journal 49(2): 379-401.
Agrawal, R. and Z. T. Fidkowski (1998). "Are thermally coupled distillation columns always
thermodynamically more efficient for ternary distillations?" Industrial & Engineering
Chemistry Research 37(8): 3444-3454.
AspenTech Aspen Properties.
Barnicki, S. D. and J. J. Siirola (2004). "Process synthesis prospective." Computers & Chemical
Engineering 28(4): 441-446.
Doherty, M. F. and M. F. Malone (2001). Conceptual Design of Distillation Systems. New York,
NY USA, McGraw-Hill.
Doherty, M. F. and M. F. Malone (2001). Conceptual design of distillation systems, McGrawHill.
Engelien, H. K. and S. Skogestad (2005). "Minimum energy diagrams for multieffect distillation
arrangements." AIChE Journal 51(6): 1714-1725.
Giridhar, A. and R. Agrawal (2010). "Synthesis of distillation configurations. II: A search
formulation for basic configurations." Computers & Chemical Engineering 34(1): 84-95.
Giridhar, A. and R. Agrawal (2010). "Synthesis of distillation configurations: I. Characteristics
of a good search space." Computers & Chemical Engineering 34(1): 73-83.
Grossmann, I. E., P. A. Aguirre, et al. (2005). "Optimal synthesis of complex distillation
columns using rigorous models." Computers & Chemical Engineering 29(6): 1203-1215.
Harwardt, A., S. Kossack, et al. (2008). Optimal column sequencing for multicomponent
mixtures. Computer Aided Chemical Engineering. B. Bertrand and J. Xavier, Elsevier.
Volume 25: 91-96.
Hilde K. Engelien, S. S. (2005). "Minimum energy diagrams for multieffect distillation
arrangements." AIChE Journal 51(6): 1714-1725.
Holland, S. T., R. Abbas, et al. (2009). "Complex Column Design by Application of Column
Profile Map Techniques: Sharp-Split Petlyuk Column Design." Industrial & Engineering
Chemistry Research 49(1): 327-349.
Holland, S. T., M. Tapp, et al. (2004). "Column profile maps. 2. Singular points and phase
diagram behaviour in ideal and nonideal systems." Industrial & Engineering Chemistry
Research 43(14): 3590-3603.
Huss, R. S. and A. W. Westerberg (1996). "Collocation Methods for Distillation Design. 1.
Model Description and Testing." Industrial & Engineering Chemistry Research 35(5):
1603-1610.
Huss, R. S. and A. W. Westerberg (1996). "Collocation Methods for Distillation Design. 2.
Applications for Distillation." Industrial & Engineering Chemistry Research 35(5): 16111623.
Kim, S. B., G. J. Ruiz, et al. (2010). "Rigorous Separation Design. 1. Multicomponent Mixtures,
Nonideal Mixtures, and Prefractionating Column Networks." Industrial & Engineering
Chemistry Research 49(14): 6499-6513.
132
Ling, H. and W. L. Luyben (2009). "New Control Structure for Divided-Wall Columns."
Industrial & Engineering Chemistry Research 48(13): 6034-6049.
Lucia, A., A. Amale, et al. (2008). "Distillation pinch points and more." Computers and
Chemical Engineering 32(6): 1350-1372.
Lucia, A. and R. Taylor (2006). "The geometry of separation boundaries: I. Basic theory and
numerical support." AIChE Journal 52(2): 582-594.
Lucia, A. and R. Taylor (2007). "The geometry of separation boundaries. II. Mathematical
formalism." AIChE Journal 53(7): 1779-1788.
Maier, R. W., J. F. Brennecke, et al. (1998). "Reliable computation of homogeneous azeotropes."
AIChE Journal 44(8): 1745-1755.
Maier, R. W., J. F. Brennecke, et al. (2000). "Reliable computation of reactive azeotropes."
Computers & Chemical Engineering 24(8): 1851-1858.
Malone, M. F. and M. F. Doherty (2000). "Reactive Distillation." Industrial & Engineering
Chemistry Research 39(11): 3953-3957.
Marquardt, W., S. Kossack, et al. (2008). "A framework for the systematic design of hybrid
separation processes." Chinese Journal of Chemical Engineering 16(3): 333-342.
Ruiz, G. J., S. Kim, et al. (2009). Design and Optimization of Energy Efficient Complex
Separation Networks. 7th International Conference on the Foundations of ComputerAided Process Design, Taylor & Francis Group.
Rump, S. M. (2009). "INTLAB." Institute of Reliable Computing http://www.ti3.tuharburg.de/.
Siirola, J. (1996). Industrial Application of Chemical Process Synthesis. Advances in chemical
engineering. J. Anderson. New York,, Academic Press.
Tapp, M., S. T. Holland, et al. (2004). "Column Profile Maps. 1. Derivation and Interpretation."
Industrial & Engineering Chemistry Research 43(2): 364-374.
Taylor, R. (2007). "(Di)Still Modeling after All These Years: A View of the State of the Art."
Industrial & Engineering Chemistry Research 46(13): 4349-4357.
Taylor, R. and R. Krishna (1993). Multicomponent Mass Transfer, Wiley.
Van Dongen, D. B. and M. F. Doherty (1985). "Design and synthesis of homogeneous azeotropic
distillations. 1. Problem formulation for a single column." Industrial & Engineering
Chemistry, Fundamentals 24(4): 454-463.
Widagdo, S. and W. D. Seider (1996). "Journal review. Azeotropic distillation." AIChE Journal
42(1): 96-130.
Wright, R. O. (1949). U.S. 2471134.
Zhang, L. and A. A. Linninger (2004). "Temperature collocation algorithm for fast and robust
distillation design." Industrial & Engineering Chemistry Research 43(12): 3163-3182.
Zhang, L. and A. A. Linninger (2006). "Towards computer-aided separation synthesis." AIChE
Journal 52(4): 1392-1409.