BACKGROUND
toward the construction of synthetic cells
one of the major goals of synthetic biology is the de novo creation of living organisms in laboratory, an ambitious and challenging goal that promises a profound impact in basic science and in biotechnology  <cit> .
the first stage towards this long-term goal is represented by the construction of quasi-cellular systems  that are liposome-based compartment containing the minimal and sufficient number of biomolecules that, by interacting with each other within such a membrane compartment, are able to emulate some of the fundamental properties of living systems.
the construction of minimal cells is carried out by following the "semi-synthetic approach" <cit> : thanks to the self-organizing behavior of lipid molecules, spherical cell-like microcompartments - called lipid vesicles  - form spontaneously in aqueous environments, encapsulating the molecular species present in solution .
typically, nucleic acids, enzymes and other biomolecules are used, in order to reconstruct minimal versions of cellular functions such as nucleic acids replication, proteins synthesis or membrane growth.
the ultimate goal of the semi-synthetic approach is the reconstruction of a cell-like system, capable of self-producing all its parts, and therefore able to grow and divide recursively.
in addition to its potential relevance for biotechnology, this approach might also shed light on the biophysics and biochemistry related to the origin of early cells on earth  <cit> .
in the past, enzymatic reactions, such as oligonucleotide synthesis  <cit>  or even pcr  <cit> , were carried inside lipid vesicle.
in  <dig> it was reported the first micro-compartmentalized polypeptide synthesis  <cit> , and soon after the first synthesis of a functional protein   <cit> .
to date, several reports deal with the fully understanding of the intra-vesicle protein synthesis, which is considered the keystone of the construction of semi-synthetic minimal cells .
in fact, the production of functional water-soluble or lipid-soluble  <cit>  proteins paves the way to the reconstruction of the minimal number of cellular functions to achieve a full self-reproducing system.
the encapsulation of solutes inside lipid vesicles
the key process underlying the construction of semi-synthetic minimal cells is the encapsulation  of biomolecules in lipid vesicles.
when such process occurs spontaneously, lipid microcompartments of several size and morphology, each containing different quantities of entrapped solutes are typically formed.
due to the stochastic nature of this encapsulation process, the concentrations of the molecular species inside nanoscale  vesicles are not always similar to the external solution.
in the case of a single solute, its intra-vesicle concentration can be lower or higher than the bulk one; in the case of multiple solutes, each species can be encapsulated with a different efficiency, creating a vast array of possibilities.
all this means that each individual vesicle behave as a unique microscopic bioreactor having a unique composition in terms of internalized components, which affects its performances in terms of effectiveness of internal reactions .
there are clear evidences that populations of spontaneously formed vesicles, irrespective of their sizes, are heterogeneous in terms of individual performances .
flow cytometry  <cit>  and direct microscopic observations  <cit>  have been used to monitor the time course of protein synthesis inside large  vesicles.
it is more difficult to perform similar studies in submicron vesicles, where it is expected, however, a higher diversity due to the enhancement of stochastic events at a smaller scale.
quite interestingly, two recent studies demonstrated that when  <dig> nm  lipid vesicles are prepared in presence of single macromolecular species, like ferritin  <cit>  or ribosomes  <cit> , the solute encapsulation does not follow the expected behavior.
it resulted, instead, in a dramatic dichotomy between many "empty" vesicles and few "super-crowded" ones .
this reveals a marked non-random behavior of the entrapment process, such that the frequency distribution of the super-crowded liposomes seems to follow the zipf-mandlebrot law that is a power-law distribution.
the volume dependency of this super-concentration effect seems to follow a power law too, resulting to be extremely marked as the vesicle diameter decreases.
these findings show that the issue of multi-solute encapsulation inside liposomes is a complex event that is far to be completely understood and that heavily affects the reactions occurring inside semi-synthetic minimal cells.
protein synthesis inside liposomes
as it has been mentioned, the technology for constructing semi-synthetic minimal cells is based on the integration of liposomes and cell-free technologies.
protein synthesis is considered one of the most relevant biological processes, and therefore it is currently studied in great detail.
most of the work is carried out by employing a completely controllable protein production system, the puresystem™   <cit> .
it contains more than  <dig> macromolecular species , and represents the minimal collection of biochemical components able to afford protein production from a coding dna sequence.
this feature is of fundamental importance considering the necessity to use known and controllable molecular components to reconstruct biological functions inside lipid micro-compartment, and moreover it is well-characterized in terms of its part, following the philosophy of synthetic biology.
thanks to the encapsulation of all ps species  inside lipid vesicles, it is possible to produce not only the gfp, but also functional enzymes as reported in several recent studies .
clearly, the vesicle encapsulation of the ps should be viewed as a complex event, that is regulated by several solute-solute and solute-membrane interactions, by stochastic  factors, and ultimately represents a realistic model of molecular self-organization at nanoscale range.
as a result of this complex pattern, a vesicle population will be characterized by functional and structural heterogeneity that need to be studied without averaging out the inter-vesicle differences.
in fact, the spreading of vesicle properties results in a kind of phenotypic diversity  <cit> .
approach
general strategy
considering the relevance of protein synthesis in the construction of semi-synthetic minimal cells, and its potential new biotechnological applications  <cit> , it is necessary to give a complete description of the protein production kinetics inside each lipid compartment.
the reasons for this need are manifold.
first, for progressing in the semi-synthetic approach to better design the construction of synthetic cells, we need to have a complete understanding of compartmentalized biochemical processes.
second, due to the difficulty of measuring the intra-vesicle composition in the case of complex mixtures, we need a tool for its evaluation, even if indirect.
third, despite the recent expansion of the field, not many studies have revealed the diversity of reaction kinetics inside vesicles and in particular in the case of the complex reaction that leads to protein synthesis.
an experimental approach uses liposomes of different size containing the dna sequence encoding for a reporter molecule  together with the ps, to detect and monitor the internal protein production over time.
in the case of relatively large vesicles  it is possible, at least in principle, to follow the course of the reaction in each vesicle thanks to light microscopy.
flow cytometry can also be used but only for following vesicles subpopulations.
large vesicles, however, due to their large size, show a minor degree in functional heterogeneity, that in turn derives from a more homogeneous internal composition.
in the case of submicron vesicles, where encapsulation anomalies are amplified  it is also possible to encapsulate the ps and produce a functional protein  <cit>  revealing, at the same time, the "conundrum" of co-encapsulation.
in this case, where multiple co-encapsulation events occur, the "anomalous entrapment" phenomenon  acts by creating nanoscale liposomes containing hundreds of different molecules in a very small volume.
in both cases , it is difficult to measure directly the exact composition of each vesicle, putting severe limitations to the full understanding of the phenomenon of multiple co-entrapment.
however, the gfp fluorescence signal emission allows the detection of the kinetics of protein production inside lipid vesicles, representing an important experimentally measurable parameter.
as we have remarked, the different rate of protein production is due to the fact that each lipid vesicle encloses a different amount of each ps species, giving rise to a unique intra-vesicle biochemical composition.
with the purpose to obtain a tool which can correctly predict the kinetics of gfp production from different intra-vesicle ps compositions, we employed a computational approach to characterize in the detail the network of biochemical reactions of the entire transcription/translation process.
computational approach
in nanoscaled biochemical environment randomness and uncertainty cannot be described simply as an additive noise factor, but they represents the fundamental forces which drive the evolution of the system  <cit> .
nevertheless, to the best of our knowledge, a detailed stochastic description of in lipo protein synthesis has discussed in a very few papers  <cit> .
to give a correct description of such small volume kinetics, the well-known gillespie stochastic simulation algorithm  <cit>  was used, implemented in a simulation software.
a model for the translation module of the ps was published by our research team  <cit> , comprising over  <dig> biochemical reactions with their kinetic coefficients; although this model is very complex, it was conceived to give a qualitative description of protein synthesis in different sized compartments, according to different entrapment models.
as every first attempt, it contained some simplifications:
1) the presence of multiple elongating ribosomes on the same rna molecule was not described;
2) protein production was conceived without taking into account aminoacids consumption;
3) the transcription process was not included;
these three problems are deeply intertwined, because the transcription and the translation processes are dynamically coupled in the ps: in the ribosomes bind to the shine-dalgarno sequence  as soon as it is available, even if the complete rna sequence has not been entirely produced yet.
this means that the translation process begins as the rna polymerase is still transcribing the rna molecule, and the two processes  cannot be treated as independent parallel processes.
moreover, different ribosomes were detected as progressing along the same rna molecule, and the same phenomenon was observed for transcription: multiple rna polymerases can simultaneously transcribe the same dna molecule.
a suitable in silico model for the ps must describe correctly all these interdependent processes.
by following the simulation strategy described in the "materials and methods" section, we obtained a new improved in silico formalization of the ps that accounts for a detailed description of the single molecular reactions; aminoacids and nucleotides must be consumed in the right quantities to ultimately produce proteins and rna.
the final goal of our approach is to investigate, by means of the obtained ps model, whether different initial concentrations of class of chemical species  can give rise to significantly different kinetics of gfp production in bulk.
then we scaled the simulated liposomes down to nano-dimensions to assess whether these differences are detectable also in such very small vesicles.
our simulation approach allows us to forecast what gfp-production kinetics will be observed in dependence of different initial concentrations.
thus our model can be used in a reverse-engineering approach where the fluorescence data experimentally recorded will be used to infer the concentration of chemical species inside liposomes.
this can reveal possible abnormalities in solutes entrapping mechanisms at nanoscale volumes.
RESULTS
ps formal model
the obtained model consists of circa  <dig> different virtual species and  <dig> reactions, among which the vast majority represents dummy species or different molecular states used to obtain sequentiality through the strategy described in the materials and methods section.
the full model is detailed in the sup.
mat., here we describe the main simplifications.
only gtp and atp are included in the model, as the two fundamental energy resources used for peptide elongation , aminoacyl-trna charging  and transcription .
the incorporation of also ctp and utp requires the explicit declaration of additional ≈  <dig> virtual species and even more biochemical reactions.
simulations experiments carried by adding ctp and utp to the presented model output results undistinguishable from those of the simplified model .
the forward motion of rna polymerases and ribosomes is dependent by the availability of dna/rna sites; this formulation does not take into account whether the available sites are on the same dna/rna strand or not.
this simplification is not affecting the consistency of the model, as the average number of elongating rna polymerases per dna strand remains the same using  <dig>   <dig> or  <dig> dna molecules.
amino acids, aminoacyl-trna synthetases, trnas, and release factors are modeled each as a single species, allowing a strong reduction of computational costs, which are notoriously high for ssa-based in silico experiments.
despite these simplifications, the portrayed model avoids the use of simplified average macroscopic measures for the formalization of initiation, elongation or termination events.
moreover, this model accounts for the presence of ordered sequential events in the elongation steps, taking into account the steric repulsions between molecules and the correct sites occupancies.
the presence of multiple elongation events, dynamical coupling between transcription and translation and sequentially ordered motion of molecules were permitted in the model formulation itself, avoiding the addition of new reactions at fixed times which can jeopardize the stochastic description of the system  <cit> .
overall consistency of the ps model
the validation of the proposed model has been obtained by comparing the rna and protein production rates with available literature data, which show a range between  <dig>  to  <dig> nt/s for the transcription process and a range between  <dig>  to  <dig> aa/s  <cit>  for translation.
extracted data for polymerases and ribosomes activities for the transcription/translation processes - <dig> nt/s and  <dig> aa/s, respectively- are consistent with experimental records.
varying the initial concentrations in large reaction volumes
different biochemical compositions were tested in silico, using the standard concentrations of species in ps reported in  <cit> , and varying the initial quantities of the  <dig> main classes of chemical species: dna , enzymes  or consumables .
they were lowered to 2/ <dig> or 1/ <dig> of their original value, for a total of  <dig> =  <dig> combinations; each combination is defined by a series of  <dig> digits, accounting for dna, enzymes or consumables concentrations; each digit contains a number which is  <dig>  or  <dig>  meaning respectively 1/ <dig>  2/ <dig> or 3/ <dig> of their normal concentration.
for example, "201" means 3/ <dig> of dna, 1/ <dig> of enzymes and 2/ <dig> of consumables.
data for protein production over time was nicely fitted  by sigmoid curves
y=gfpmax1+e-x-x <dig> gfpmaxb
with parameters:
gfpmax = maximum value of protein produced
b = maximum slope of the curve
x <dig> gfpmax=time value for y=gfpmax/ <dig>
for each initial biochemical combination we extracted the gfpmax and b parameters, to compare the total gfp yield and its production rate for the different ps compositions .
general dependencies by the initial conditions
high dna concentrations accelerates the overall protein production; however, the overall protein yield  diminishes as the dna amount increases; in fact, although rapid, protein production with high dna concentrations stops at lower time values .
simulations carried with a lower amount of enzymes concentrations resulted in a strong decrease of protein yield: a change from 3/ <dig> to 2/ <dig> in enzymes concentration determined a reduction in protein yield to circa 40% of the total; an additional reduction in enzymes concentration to a 1/ <dig> of the original value led to a 5% of produced protein compared to normal conditions.
in addition, very slow kinetics  are observed as the enzymes concentration decreases.
input files with 1/ <dig> of initial consumables amount yielded a maximum of only  <dig> protein molecules , revealing the role of energy resources as the fundamental factor which most affects the overall gfp production.
the highest protein production is afforded when dna is low and enzyme and consumables are present in maximum quantity , reaching a total of approximately  <dig> gfp molecules .
trying to establish the exact role of dna concentration in the protein production kinetics we performed stochastic simulations experiments using different lower dna initial amounts.
as reported before, low dna concentrations determine a slow kinetic of protein production: simulation performed using very low dna amounts resulted in unfeasibly slow gfp productions.
furthermore, protein production encompasses several hours of time, and self-inactivating phenomena, which probably involve ribosomes inactivation, were experimentally observed after approximately  <dig> hours from the beginning of the experiment  <cit> ;
thus, ad hoc simulations were performed to account for a self-inactivating reaction for ribosomes after circa  <dig> hours.
data was fitted by sigmoid curves and the extracted parameters were compared as discussed before.
final protein yields are shown for different initial dna amounts.
as dna concentration decreases, the transcription process slows down, and gfp production kinetics are slower, but the total protein yield is increased to a total of over  <dig>  gfp molecules  when  = 15nm, a result which is comparable with experimental measures in gfp-expressing giant lipid vesicles of equal volume  <cit> .
the non-linear dependence of total protein yield on the initial dna concentration is the result of the balancing of two opposite phenomena: at very low initial dna concentrations, the protein production goes slow, but almost all the chemical energy contained in the liposome can be used to produce the protein.
on the other hand, when initial dna concentration is high, the protein production is fast, but the energy consumption is even faster, making it impossible to obtain a larger protein yield .
in fact, these results about the behavior of the ps  showed how the competition between the ribosomes inactivation process and protein production speed determines a critical dna concentration value, which delineates an optimal distribution of energy resources between the different processes .
energetic assessment of the ps
when using lower initial dna concentrations  the competition effect between the transcription and translation processes for energy resources becomes a clearer phenomenon: figure  <dig> illustrates the primal role of dna in determining the overall gfp yield, and how scaling the initial amount of  results in a different energy consumption between the biochemical processes present in the ps reaction network.
when dna concentration is high, transcription produces immediately many ribosomal binding sites and the translation initiation  is more likely to occur, resulting in a more rapid protein production kinetic compared with inputs containing lower dna concentrations, but at the same time, the transcription process consumes a large amount of nucleotides during the rna elongation process.
translation factors  use gtp as energy donor, thus the lack of gtp molecules causes the translation process, and subsequently protein production, to stop.
ndk  can provide new gtp molecules from the atp pool, but the additional gtp is consumed by transcription or in the intermediate translation steps, unlikely resulting in the formation of a significant number of new complete proteins.
this simulation study allows us to assess two results.
firstly, the small changes in the initial concentrations of the ps chemical classes significantly influence the final gfp production and kinetics .
secondly, the initial dna concentration can be deduced from the global gfp production, thus using the proposed model for a reverse engineering approach.
at nanoscale range
the different protein yields for the different ps combinations were tested also in small volumes ; parameters were extracted after data fitting using the same procedure discussed before.
the obtained results  showed a similar distribution of the gfp production efficiencies, reaching a maximum of approximately  <dig> proteins  for the "022" sample, as seen also for simulations in higher volumes.
many combinations resulted in negligible gfp production, indicating the importance in lower volume of an optimal internal biochemical composition.
the general trend for protein yield/kinetics and its dependencies by dna/enzymes/consumables concentrations is again predicted in small volumes, showing the behavior of the ps when no anomalous entrapment phenomena are supposed.
in all the  <dig> different initial conditions, we can observe the relative low importance of the dna quantity, probably because it is anyhow relatively high with respect to the initial concentration of enzymes and consumables.
this situation is a direct consequence of the high discretization level: to have a fully functional liposome, at least one complete dna molecule should be entrapped.
this may alter the exact proportionality with the initial concentration in bulk.
we finally observe that the combination "x22"  shows the maximum estimated protein production.
to suggest an experimental test, we recommend to compare the "200" and the "022" samples, in order to have the maximum possible resolution for the differences in gfp production, between the two initial conditions.
concerning the kinetics studies, we confirm that the two optimal experimental candidates are the "200" and "022" samples, but we observe that the kinetic parameter is critically dependent on the initial enzymes concentration, that, at these small volumes, can show very high random fluctuations .
we remark that also in small volumes, the small changes in initial concentrations lead to significantly different final statuses, thus supporting the proposed experimental design.
CONCLUSIONS
with this work we have conquered a methodological goal, finding a way to describe sequential processes in a concurrent-only stochastic simulation environment.
thanks to this we have given the most detailed description of the ps ever proposed.
we have then explored different initial concentration and forecasted the kinetics of protein synthesis for large and small volume liposomes.
within the small volumes, we isolated extreme final values of gfp-production in order to select candidate of experimentally distinguishable cases.
we present our platform as a reverse-engineering tool to be used to analyze future experimental data.
materials and methods
qdc and its input language
the formal specification of the ps reaction network was assessed using a stochastic simulation software previously built by our team, qdc , which is based on the direct method version of the gillespie stochastic simulation algorithm  <cit> .
qdc's input language is simple and intuitive, and easily recalls the standard notation for biochemical reactions.
qdc has been designed to simulate experiments performed on metabolic networks, where the operator can exert three different control action on the metabolic network: 1) add or remove molecules at a given time; 2) change the propensity of a reaction at a given time; 3) simulate all-or-nothing reactions , where the chemical reaction is executed immediately after the stoichiometry of the left part of the equation is satisfied .
immediate reactions are not standard biochemical events that occur accordingly with a kinetic law, but they rather represent logical statements, which were introduced to allow the description of complex condition with possibly many chemical species interacting.
these features of qdc allowed us to build a detailed model of the ps reaction network.
biochemical reactions with their relative constants from the previous model  <cit>  were updated, when possible, using data from literature and the bionumbers  <cit>  database.
the transcription reactions were modeled using data from published material  <cit> ; the translational core model was updated using a more detailed kinetic description  <cit> .
in silico puresystem model
the first result obtained in this work has been a detailed formal description of the ps: every molecular component was explicitly inserted in the model together with the known biochemical reactions between the different reagents, as described in the previous section.
the same model simplifications present in the previous work  <cit>  were introduced to create a suitable transcription/translation formulation.
however, as pointed in the previous section, a formal description of interdependent elongation events transcends the standard ssa specifications, which depicts the system evolution according to propensities rules only proportional to the number of reacting molecules, without accounting for the presence of sequential processes or topological constraints .
to coherently describe the proper distance between elongating molecules and, subsequently, a correct occupancy of rna and dna sequences, all the species involved in the elongation events  were split into multiple entities representing different molecular states .
the divisions in multiple molecular states allows a coherent organization of site occupancies and dynamic motion, but the transition events describing the forward motion of ribosomes and polymerases are two processes which act in a sequential fashion.
the solution to this problem was assessed using logical statements written in form of reactions present in the qdc syntax, the aforementioned immediate reactions, which were massively used in this new model to regulate the forward motion of polymerases and ribosomes, accounting for correct spacing and dynamical coupling.
modeling strategy
the dna sequence encoding for gfp was divided, according its length, into different multiple species, each representing a  <dig> bp sequence; the polymerization process is divided into different reactions, describing a second-order reaction for nucleotide binding, and a first-order reaction for nucleotides incorporation which returns the polymerase molecule  and a "dummy" product, that allows to track the number of nucleotides incorporated in the rna molecule ; the immediate reactions determine the transition to the next step, ensuring the following conditions: a) an adjacent dna site is available, b) a correct number of nucleotides has been added to the rna sequence, c) the corresponding rna sequence is produced, d) the previously occupied dna site is released.
here it is an example for transcription process at the fourth step :
<dig>  t7elgt <dig> + gtp > t7pregel4
<dig>  t7pregel <dig> > t7elat <dig> + pi + gtr4
<dig>  t7elat <dig> + atp > t7preael4
<dig>  t7preael <dig> >t7elgt <dig> + pi + atr4
-,  <dig> gtr <dig> +  <dig> atr <dig> + t7elgt <dig> + dna <dig> > t7elgt <dig> + rna <dig> + dna3
this reaction "box" was duplicated different times ensuring the correct succession of molecular states; when only one dna molecule is available, the polymerases advance one by one, separated by at least one dna site between each other .
the same strategy was used to describe the translation reactions.
<dig>  er <dig> + efargtp > expre3
<dig>  expre <dig> > er <dig> + efargtp
<dig>  expre <dig> > ex3
<dig> , ex <dig> > expre3
<dig>  ex <dig> > ex3gdp + pi
<dig>  ex3gdp > era3gdp
<dig>  era3gdp > era <dig> + eftugdp
<dig>  era <dig> > era3trna
<dig>  era3trna + efggtp > exbpre3
<dig>  exbpre <dig> > era3trna + efggtp
<dig>  exbpre <dig> > exb <dig> + efg + gdp +pi
<dig>  exb <dig> >er <dig> + trna + transl3
-,  <dig> transl <dig> + rna <dig> + pept <dig> + er <dig> > er <dig> + pept <dig> + rna2
an elongating ribosome  binds the complex which carry the aminoacid , after which moves to the next codon, aided by the elongation factor efg charged with gtp ; this translocation reaction yield ad additional product which is used to regulate the progression to the next state.
after a fixed number of translocation steps  an immediate reaction occurs in a similar fashion as seen for transcription: 1) the correct amount of aminoacids are incorporated, and thus consumed; 2) the next free rna site is occupied and 3) the previous one is therefore liberated; 4) an entity named pept is also produced, allowing to calculate the length of the peptide sequence produced so far: for example, if  <dig> species named pept <dig> are present in a certain time of the simulation, this means that there are  <dig> peptides, still bound to the ribosomes, with a length spanning from  <dig> ×  <dig> =  <dig> to - <dig> =  <dig> aminoacids.
additional file  <dig> contains a legend of all the chemical species declared in our in silico ps, while additional file  <dig> contains the complete input file for simulating by qdc the ps entrapped in a 10- <dig> l liposome.
list of abbreviations used
ps: puresystem™; popc: 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphatidylcholine; ndk: nucleotide diphosphate kinase; gfp: green fluorescent protein; atp: adenosine triphosphate; gtp: guanosine triphosphate.
competing interests
the authors declare that they have no competing interests.
supplementary material
additional file 1
a text file with the legend of all the names used for chemical species in the input file.
click here for file
additional file 2
the input file for a liposome of 10- <dig> l; this is readable by qdc simulator.
click here for file
additional file 3
a ms-excel file containing all the results of the simulations performed for large-volume liposomes.
click here for file
additional file 4
a ms-excel file containing all the results of the simulations performed for very low initial dna concentrations.
click here for file
additional file 5
a ms-excel file containing all the results of the simulations performed for small-volume liposomes.
click here for file
acknowledgements
the authors thank claudio felicioli for his valuable collaboration in managing the code of the stochastic simulator qdc.
part of this work was developed within the european cost systems chemistry action .
declarations
the authors declare that this paper is part of the supplement dedicated to bits <dig> conference.
the funds covering the publication fee have been collected from research funds assigned by authors' departments and the miur project prin 2008fy7rj <dig>
this article has been published as part of bmc bioinformatics volume  <dig> supplement  <dig>  2013: italian society of bioinformatics : annual meeting  <dig>  the full contents of the supplement are available online at http://www.biomedcentral.com/bmcbioinformatics/supplements/14/s7

 toward the construction of synthetic cells one of the major goals of synthetic biology is the de novo creation of living organisms in laboratory, an ambitious and challenging goal that promises a profound impact in basic science and in biotechnology. the first stage towards this long-term goal is represented by the construction of quasi-cellular systems that are liposome-based compartment containing the minimal and sufficient number of biomolecules that, by interacting with each other within such a membrane compartment, are able to emulate some of the fundamental properties of living systems. the construction of minimal cells is carried out by following the "semi-synthetic approach": thanks to the self-organizing behavior of lipid molecules, spherical cell-like microcompartments - called lipid vesicles - form spontaneously in aqueous environments, encapsulating the molecular species present in solution. typically, nucleic acids, enzymes and other biomolecules are used, in order to reconstruct minimal versions of cellular functions such as nucleic acids replication, proteins synthesis or membrane growth. the ultimate goal of the semi-synthetic approach is the reconstruction of a cell-like system, capable of self-producing all its parts, and therefore able to grow and divide recursively. in addition to its potential relevance for biotechnology, this approach might also shed light on the biophysics and biochemistry related to the origin of early cells on earth. in the past, enzymatic reactions, such as oligonucleotide synthesis or even pcr, were carried inside lipid vesicle. in it was reported the first micro-compartmentalized polypeptide synthesis, and soon after the first synthesis of a functional protein. to date, several reports deal with the fully understanding of the intra-vesicle protein synthesis, which is considered the keystone of the construction of semi-synthetic minimal cells. in fact, the production of functional water-soluble or lipid-soluble proteins paves the way to the reconstruction of the minimal number of cellular functions to achieve a full self-reproducing system. the encapsulation of solutes inside lipid vesicles the key process underlying the construction of semi-synthetic minimal cells is the encapsulation of biomolecules in lipid vesicles. when such process occurs spontaneously, lipid microcompartments of several size and morphology, each containing different quantities of entrapped solutes are typically formed. due to the stochastic nature of this encapsulation process, the concentrations of the molecular species inside nanoscale vesicles are not always similar to the external solution. in the case of a single solute, its intra-vesicle concentration can be lower or higher than the bulk one; in the case of multiple solutes, each species can be encapsulated with a different efficiency, creating a vast array of possibilities. all this means that each individual vesicle behave as a unique microscopic bioreactor having a unique composition in terms of internalized components, which affects its performances in terms of effectiveness of internal reactions. there are clear evidences that populations of spontaneously formed vesicles, irrespective of their sizes, are heterogeneous in terms of individual performances. flow cytometry and direct microscopic observations have been used to monitor the time course of protein synthesis inside large vesicles. it is more difficult to perform similar studies in submicron vesicles, where it is expected, however, a higher diversity due to the enhancement of stochastic events at a smaller scale. quite interestingly, two recent studies demonstrated that when nm lipid vesicles are prepared in presence of single macromolecular species, like ferritin or ribosomes, the solute encapsulation does not follow the expected behavior. it resulted, instead, in a dramatic dichotomy between many "empty" vesicles and few "super-crowded" ones. this reveals a marked non-random behavior of the entrapment process, such that the frequency distribution of the super-crowded liposomes seems to follow the zipf-mandlebrot law that is a power-law distribution. the volume dependency of this super-concentration effect seems to follow a power law too, resulting to be extremely marked as the vesicle diameter decreases. these findings show that the issue of multi-solute encapsulation inside liposomes is a complex event that is far to be completely understood and that heavily affects the reactions occurring inside semi-synthetic minimal cells. protein synthesis inside liposomes as it has been mentioned, the technology for constructing semi-synthetic minimal cells is based on the integration of liposomes and cell-free technologies. protein synthesis is considered one of the most relevant biological processes, and therefore it is currently studied in great detail. most of the work is carried out by employing a completely controllable protein production system, the puresystem™. it contains more than macromolecular species, and represents the minimal collection of biochemical components able to afford protein production from a coding dna sequence. this feature is of fundamental importance considering the necessity to use known and controllable molecular components to reconstruct biological functions inside lipid micro-compartment, and moreover it is well-characterized in terms of its part, following the philosophy of synthetic biology. thanks to the encapsulation of all ps species inside lipid vesicles, it is possible to produce not only the gfp, but also functional enzymes as reported in several recent studies. clearly, the vesicle encapsulation of the ps should be viewed as a complex event, that is regulated by several solute-solute and solute-membrane interactions, by stochastic factors, and ultimately represents a realistic model of molecular self-organization at nanoscale range. as a result of this complex pattern, a vesicle population will be characterized by functional and structural heterogeneity that need to be studied without averaging out the inter-vesicle differences. in fact, the spreading of vesicle properties results in a kind of phenotypic diversity. approach general strategy considering the relevance of protein synthesis in the construction of semi-synthetic minimal cells, and its potential new biotechnological applications, it is necessary to give a complete description of the protein production kinetics inside each lipid compartment. the reasons for this need are manifold. first, for progressing in the semi-synthetic approach to better design the construction of synthetic cells, we need to have a complete understanding of compartmentalized biochemical processes. second, due to the difficulty of measuring the intra-vesicle composition in the case of complex mixtures, we need a tool for its evaluation, even if indirect. third, despite the recent expansion of the field, not many studies have revealed the diversity of reaction kinetics inside vesicles and in particular in the case of the complex reaction that leads to protein synthesis. an experimental approach uses liposomes of different size containing the dna sequence encoding for a reporter molecule together with the ps, to detect and monitor the internal protein production over time. in the case of relatively large vesicles it is possible, at least in principle, to follow the course of the reaction in each vesicle thanks to light microscopy. flow cytometry can also be used but only for following vesicles subpopulations. large vesicles, however, due to their large size, show a minor degree in functional heterogeneity, that in turn derives from a more homogeneous internal composition. in the case of submicron vesicles, where encapsulation anomalies are amplified it is also possible to encapsulate the ps and produce a functional protein revealing, at the same time, the "conundrum" of co-encapsulation. in this case, where multiple co-encapsulation events occur, the "anomalous entrapment" phenomenon acts by creating nanoscale liposomes containing hundreds of different molecules in a very small volume. in both cases, it is difficult to measure directly the exact composition of each vesicle, putting severe limitations to the full understanding of the phenomenon of multiple co-entrapment. however, the gfp fluorescence signal emission allows the detection of the kinetics of protein production inside lipid vesicles, representing an important experimentally measurable parameter. as we have remarked, the different rate of protein production is due to the fact that each lipid vesicle encloses a different amount of each ps species, giving rise to a unique intra-vesicle biochemical composition. with the purpose to obtain a tool which can correctly predict the kinetics of gfp production from different intra-vesicle ps compositions, we employed a computational approach to characterize in the detail the network of biochemical reactions of the entire transcription/translation process. computational approach in nanoscaled biochemical environment randomness and uncertainty cannot be described simply as an additive noise factor, but they represents the fundamental forces which drive the evolution of the system. nevertheless, to the best of our knowledge, a detailed stochastic description of in lipo protein synthesis has discussed in a very few papers. to give a correct description of such small volume kinetics, the well-known gillespie stochastic simulation algorithm was used, implemented in a simulation software. a model for the translation module of the ps was published by our research team, comprising over biochemical reactions with their kinetic coefficients; although this model is very complex, it was conceived to give a qualitative description of protein synthesis in different sized compartments, according to different entrapment models. as every first attempt, it contained some simplifications: 1) the presence of multiple elongating ribosomes on the same rna molecule was not described; 2) protein production was conceived without taking into account aminoacids consumption; 3) the transcription process was not included; these three problems are deeply intertwined, because the transcription and the translation processes are dynamically coupled in the ps: in the ribosomes bind to the shine-dalgarno sequence as soon as it is available, even if the complete rna sequence has not been entirely produced yet. this means that the translation process begins as the rna polymerase is still transcribing the rna molecule, and the two processes cannot be treated as independent parallel processes. moreover, different ribosomes were detected as progressing along the same rna molecule, and the same phenomenon was observed for transcription: multiple rna polymerases can simultaneously transcribe the same dna molecule. a suitable in silico model for the ps must describe correctly all these interdependent processes. by following the simulation strategy described in the "materials and methods" section, we obtained a new improved in silico formalization of the ps that accounts for a detailed description of the single molecular reactions; aminoacids and nucleotides must be consumed in the right quantities to ultimately produce proteins and rna. the final goal of our approach is to investigate, by means of the obtained ps model, whether different initial concentrations of class of chemical species can give rise to significantly different kinetics of gfp production in bulk. then we scaled the simulated liposomes down to nano-dimensions to assess whether these differences are detectable also in such very small vesicles. our simulation approach allows us to forecast what gfp-production kinetics will be observed in dependence of different initial concentrations. thus our model can be used in a reverse-engineering approach where the fluorescence data experimentally recorded will be used to infer the concentration of chemical species inside liposomes. this can reveal possible abnormalities in solutes entrapping mechanisms at nanoscale volumes.  ps formal model the obtained model consists of circa different virtual species and reactions, among which the vast majority represents dummy species or different molecular states used to obtain sequentiality through the strategy described in the materials and methods section. the full model is detailed in the sup. mat., here we describe the main simplifications. only gtp and atp are included in the model, as the two fundamental energy resources used for peptide elongation, aminoacyl-trna charging and transcription. the incorporation of also ctp and utp requires the explicit declaration of additional ≈ virtual species and even more biochemical reactions. simulations experiments carried by adding ctp and utp to the presented model output results undistinguishable from those of the simplified model. the forward motion of rna polymerases and ribosomes is dependent by the availability of dna/rna sites; this formulation does not take into account whether the available sites are on the same dna/rna strand or not. this simplification is not affecting the consistency of the model, as the average number of elongating rna polymerases per dna strand remains the same using or dna molecules. amino acids, aminoacyl-trna synthetases, trnas, and release factors are modeled each as a single species, allowing a strong reduction of computational costs, which are notoriously high for ssa-based in silico experiments. despite these simplifications, the portrayed model avoids the use of simplified average macroscopic measures for the formalization of initiation, elongation or termination events. moreover, this model accounts for the presence of ordered sequential events in the elongation steps, taking into account the steric repulsions between molecules and the correct sites occupancies. the presence of multiple elongation events, dynamical coupling between transcription and translation and sequentially ordered motion of molecules were permitted in the model formulation itself, avoiding the addition of new reactions at fixed times which can jeopardize the stochastic description of the system. overall consistency of the ps model the validation of the proposed model has been obtained by comparing the rna and protein production rates with available literature data, which show a range between to nt/s for the transcription process and a range between to aa/s for translation. extracted data for polymerases and ribosomes activities for the transcription/translation processes - nt/s and aa/s, respectively- are consistent with experimental records. varying the initial concentrations in large reaction volumes different biochemical compositions were tested in silico, using the standard concentrations of species in ps reported in, and varying the initial quantities of the main classes of chemical species: dna, enzymes or consumables. they were lowered to 2/ or 1/ of their original value, for a total of = combinations; each combination is defined by a series of digits, accounting for dna, enzymes or consumables concentrations; each digit contains a number which is or meaning respectively 1/ 2/ or 3/ of their normal concentration. for example, "201" means 3/ of dna, 1/ of enzymes and 2/ of consumables. data for protein production over time was nicely fitted by sigmoid curves y=gfpmax1+e-x-x gfpmaxb with parameters: gfpmax = maximum value of protein produced b = maximum slope of the curve x gfpmax=time value for y=gfpmax/ for each initial biochemical combination we extracted the gfpmax and b parameters, to compare the total gfp yield and its production rate for the different ps compositions. general dependencies by the initial conditions high dna concentrations accelerates the overall protein production; however, the overall protein yield diminishes as the dna amount increases; in fact, although rapid, protein production with high dna concentrations stops at lower time values. simulations carried with a lower amount of enzymes concentrations resulted in a strong decrease of protein yield: a change from 3/ to 2/ in enzymes concentration determined a reduction in protein yield to circa 40% of the total; an additional reduction in enzymes concentration to a 1/ of the original value led to a 5% of produced protein compared to normal conditions. in addition, very slow kinetics are observed as the enzymes concentration decreases. input files with 1/ of initial consumables amount yielded a maximum of only protein molecules, revealing the role of energy resources as the fundamental factor which most affects the overall gfp production. the highest protein production is afforded when dna is low and enzyme and consumables are present in maximum quantity, reaching a total of approximately gfp molecules. trying to establish the exact role of dna concentration in the protein production kinetics we performed stochastic simulations experiments using different lower dna initial amounts. as reported before, low dna concentrations determine a slow kinetic of protein production: simulation performed using very low dna amounts resulted in unfeasibly slow gfp productions. furthermore, protein production encompasses several hours of time, and self-inactivating phenomena, which probably involve ribosomes inactivation, were experimentally observed after approximately hours from the beginning of the experiment; thus, ad hoc simulations were performed to account for a self-inactivating reaction for ribosomes after circa hours. data was fitted by sigmoid curves and the extracted parameters were compared as discussed before. final protein yields are shown for different initial dna amounts. as dna concentration decreases, the transcription process slows down, and gfp production kinetics are slower, but the total protein yield is increased to a total of over gfp molecules when = 15nm, a result which is comparable with experimental measures in gfp-expressing giant lipid vesicles of equal volume. the non-linear dependence of total protein yield on the initial dna concentration is the result of the balancing of two opposite phenomena: at very low initial dna concentrations, the protein production goes slow, but almost all the chemical energy contained in the liposome can be used to produce the protein. on the other hand, when initial dna concentration is high, the protein production is fast, but the energy consumption is even faster, making it impossible to obtain a larger protein yield. in fact, these results about the behavior of the ps showed how the competition between the ribosomes inactivation process and protein production speed determines a critical dna concentration value, which delineates an optimal distribution of energy resources between the different processes. energetic assessment of the ps when using lower initial dna concentrations the competition effect between the transcription and translation processes for energy resources becomes a clearer phenomenon: figure illustrates the primal role of dna in determining the overall gfp yield, and how scaling the initial amount of results in a different energy consumption between the biochemical processes present in the ps reaction network. when dna concentration is high, transcription produces immediately many ribosomal binding sites and the translation initiation is more likely to occur, resulting in a more rapid protein production kinetic compared with inputs containing lower dna concentrations, but at the same time, the transcription process consumes a large amount of nucleotides during the rna elongation process. translation factors use gtp as energy donor, thus the lack of gtp molecules causes the translation process, and subsequently protein production, to stop. ndk can provide new gtp molecules from the atp pool, but the additional gtp is consumed by transcription or in the intermediate translation steps, unlikely resulting in the formation of a significant number of new complete proteins. this simulation study allows us to assess two results. firstly, the small changes in the initial concentrations of the ps chemical classes significantly influence the final gfp production and kinetics. secondly, the initial dna concentration can be deduced from the global gfp production, thus using the proposed model for a reverse engineering approach. at nanoscale range the different protein yields for the different ps combinations were tested also in small volumes; parameters were extracted after data fitting using the same procedure discussed before. the obtained results showed a similar distribution of the gfp production efficiencies, reaching a maximum of approximately proteins for the "022" sample, as seen also for simulations in higher volumes. many combinations resulted in negligible gfp production, indicating the importance in lower volume of an optimal internal biochemical composition. the general trend for protein yield/kinetics and its dependencies by dna/enzymes/consumables concentrations is again predicted in small volumes, showing the behavior of the ps when no anomalous entrapment phenomena are supposed. in all the different initial conditions, we can observe the relative low importance of the dna quantity, probably because it is anyhow relatively high with respect to the initial concentration of enzymes and consumables. this situation is a direct consequence of the high discretization level: to have a fully functional liposome, at least one complete dna molecule should be entrapped. this may alter the exact proportionality with the initial concentration in bulk. we finally observe that the combination "x22" shows the maximum estimated protein production. to suggest an experimental test, we recommend to compare the "200" and the "022" samples, in order to have the maximum possible resolution for the differences in gfp production, between the two initial conditions. concerning the kinetics studies, we confirm that the two optimal experimental candidates are the "200" and "022" samples, but we observe that the kinetic parameter is critically dependent on the initial enzymes concentration, that, at these small volumes, can show very high random fluctuations. we remark that also in small volumes, the small changes in initial concentrations lead to significantly different final statuses, thus supporting the proposed experimental design.  with this work we have conquered a methodological goal, finding a way to describe sequential processes in a concurrent-only stochastic simulation environment. thanks to this we have given the most detailed description of the ps ever proposed. we have then explored different initial concentration and forecasted the kinetics of protein synthesis for large and small volume liposomes. within the small volumes, we isolated extreme final values of gfp-production in order to select candidate of experimentally distinguishable cases. we present our platform as a reverse-engineering tool to be used to analyze future experimental data. materials and methods qdc and its input language the formal specification of the ps reaction network was assessed using a stochastic simulation software previously built by our team, qdc, which is based on the direct method version of the gillespie stochastic simulation algorithm. qdc's input language is simple and intuitive, and easily recalls the standard notation for biochemical reactions. qdc has been designed to simulate experiments performed on metabolic networks, where the operator can exert three different control action on the metabolic network: 1) add or remove molecules at a given time; 2) change the propensity of a reaction at a given time; 3) simulate all-or-nothing reactions, where the chemical reaction is executed immediately after the stoichiometry of the left part of the equation is satisfied. immediate reactions are not standard biochemical events that occur accordingly with a kinetic law, but they rather represent logical statements, which were introduced to allow the description of complex condition with possibly many chemical species interacting. these features of qdc allowed us to build a detailed model of the ps reaction network. biochemical reactions with their relative constants from the previous model were updated, when possible, using data from literature and the bionumbers database. the transcription reactions were modeled using data from published material; the translational core model was updated using a more detailed kinetic description. in silico puresystem model the first result obtained in this work has been a detailed formal description of the ps: every molecular component was explicitly inserted in the model together with the known biochemical reactions between the different reagents, as described in the previous section. the same model simplifications present in the previous work were introduced to create a suitable transcription/translation formulation. however, as pointed in the previous section, a formal description of interdependent elongation events transcends the standard ssa specifications, which depicts the system evolution according to propensities rules only proportional to the number of reacting molecules, without accounting for the presence of sequential processes or topological constraints. to coherently describe the proper distance between elongating molecules and, subsequently, a correct occupancy of rna and dna sequences, all the species involved in the elongation events were split into multiple entities representing different molecular states. the divisions in multiple molecular states allows a coherent organization of site occupancies and dynamic motion, but the transition events describing the forward motion of ribosomes and polymerases are two processes which act in a sequential fashion. the solution to this problem was assessed using logical statements written in form of reactions present in the qdc syntax, the aforementioned immediate reactions, which were massively used in this new model to regulate the forward motion of polymerases and ribosomes, accounting for correct spacing and dynamical coupling. modeling strategy the dna sequence encoding for gfp was divided, according its length, into different multiple species, each representing a bp sequence; the polymerization process is divided into different reactions, describing a second-order reaction for nucleotide binding, and a first-order reaction for nucleotides incorporation which returns the polymerase molecule and a "dummy" product, that allows to track the number of nucleotides incorporated in the rna molecule; the immediate reactions determine the transition to the next step, ensuring the following conditions: a) an adjacent dna site is available, b) a correct number of nucleotides has been added to the rna sequence, c) the corresponding rna sequence is produced, d) the previously occupied dna site is released. here it is an example for transcription process at the fourth step: t7elgt + gtp > t7pregel4 t7pregel > t7elat + pi + gtr4 t7elat + atp > t7preael4 t7preael >t7elgt + pi + atr4 -, gtr + atr + t7elgt + dna > t7elgt + rna + dna3 this reaction "box" was duplicated different times ensuring the correct succession of molecular states; when only one dna molecule is available, the polymerases advance one by one, separated by at least one dna site between each other. the same strategy was used to describe the translation reactions. er + efargtp > expre3 expre > er + efargtp expre > ex3, ex > expre3 ex > ex3gdp + pi ex3gdp > era3gdp era3gdp > era + eftugdp era > era3trna era3trna + efggtp > exbpre3 exbpre > era3trna + efggtp exbpre > exb + efg + gdp +pi exb >er + trna + transl3 -, transl + rna + pept + er > er + pept + rna2 an elongating ribosome binds the complex which carry the aminoacid, after which moves to the next codon, aided by the elongation factor efg charged with gtp; this translocation reaction yield ad additional product which is used to regulate the progression to the next state. after a fixed number of translocation steps an immediate reaction occurs in a similar fashion as seen for transcription: 1) the correct amount of aminoacids are incorporated, and thus consumed; 2) the next free rna site is occupied and 3) the previous one is therefore liberated; 4) an entity named pept is also produced, allowing to calculate the length of the peptide sequence produced so far: for example, if species named pept are present in a certain time of the simulation, this means that there are peptides, still bound to the ribosomes, with a length spanning from × = to - = aminoacids. additional file contains a legend of all the chemical species declared in our in silico ps, while additional file contains the complete input file for simulating by qdc the ps entrapped in a 10- l liposome. list of abbreviations used ps: puresystem™; popc: 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphatidylcholine; ndk: nucleotide diphosphate kinase; gfp: green fluorescent protein; atp: adenosine triphosphate; gtp: guanosine triphosphate. competing interests the authors declare that they have no competing interests. supplementary material additional file 1 a text file with the legend of all the names used for chemical species in the input file. click here for file additional file 2 the input file for a liposome of 10- l; this is readable by qdc simulator. click here for file additional file 3 a ms-excel file containing all the results of the simulations performed for large-volume liposomes. click here for file additional file 4 a ms-excel file containing all the results of the simulations performed for very low initial dna concentrations. click here for file additional file 5 a ms-excel file containing all the results of the simulations performed for small-volume liposomes. click here for file acknowledgements the authors thank claudio felicioli for his valuable collaboration in managing the code of the stochastic simulator qdc. part of this work was developed within the european cost systems chemistry action. declarations the authors declare that this paper is part of the supplement dedicated to bits conference. the funds covering the publication fee have been collected from research funds assigned by authors' departments and the miur project prin 2008fy7rj this article has been published as part of bmc bioinformatics volume supplement 2013: italian society of bioinformatics: annual meeting the full contents of the supplement are available online at http://www.biomedcentral.com/bmcbioinformatics/supplements/14/s7
qdc, a stochastic simulation software based on the well-known gillespie's ssa algorithm, was used. quasi-cellular systems able to produce proteins directly from dna can be obtained by encapsulating the cell-free transcription/translation system puresystem™ in liposomes. by varying the initial concentrations of the three main classes of molecules involved in the ps, we were able to stochastically simulate the time-course of gfp-production. then we extended this study for small-volume liposomes, where it is more complex to infer the intra-vesicle composition, due to the expected anomalous entrapment phenomena. thanks to our results, an experimental approach is now possible, aimed at recording the gfp production kinetics in very small micro-emulsion droplets or liposomes, and inferring, by using the simulation as a reverse-engineering procedure, the internal solutes distribution, and shed light on the still unknown forces driving the entrapment phenomenon. thanks to suitable features specific of qdc, we successfully formalized the dynamical coupling between the transcription and translation processes that occurs in the real ps, thus bypassing the concurrent-only environment of gillespie's algorithm. 2- may ninth annual meeting of the italian society of bioinformatics catania, sicily stochastic simulation is a valuable tool in the study of biochemical reaction at nanoscale range. the entrapment of solutes in small-volume liposomes is a fundamental open problem.
 the wet-lab synthesis of the simplest forms of life is a challenging aspect in modern synthetic biology. quasi-cellular systems able to produce proteins directly from dna can be obtained by encapsulating the cell-free transcription/translation system puresystem™ in liposomes. it is possible to detect the intra-vesicle protein production using dna encoding for gfp and monitoring the fluorescence emission over time. the entrapment of solutes in small-volume liposomes is a fundamental open problem. stochastic simulation is a valuable tool in the study of biochemical reaction at nanoscale range. qdc, a stochastic simulation software based on the well-known gillespie's ssa algorithm, was used. a suitable model formally describing the ps reactions network was developed, to predict, from inner species concentrations, the resulting fluorescence signal.  thanks to suitable features specific of qdc, we successfully formalized the dynamical coupling between the transcription and translation processes that occurs in the real ps, thus bypassing the concurrent-only environment of gillespie's algorithm. simulations were firstly performed for large liposomes entrapping the ps to synthetize gfp. by varying the initial concentrations of the three main classes of molecules involved in the ps, we were able to stochastically simulate the time-course of gfp-production. the sigmoid fit of the gfp-production curves allowed us to extract three quantitative parameters which are significantly dependent on the various initial states. then we extended this study for small-volume liposomes, where it is more complex to infer the intra-vesicle composition, due to the expected anomalous entrapment phenomena. we identified almost two extreme states that are forecasted to give rise to significantly different experimental observables.  the present work is the first one describing in the detail the stochastic behavior of the ps. thanks to our results, an experimental approach is now possible, aimed at recording the gfp production kinetics in very small micro-emulsion droplets or liposomes, and inferring, by using the simulation as a reverse-engineering procedure, the internal solutes distribution, and shed light on the still unknown forces driving the entrapment phenomenon. 2- may ninth annual meeting of the italian society of bioinformatics catania, sicily 