C Jayyosi

Karine Bruyère-Garnier

M Coret

Geometry of an inflated membrane in elliptic bulge tests: evaluation of an ellipsoidal shape approximation by stereoscopic digital image correlation measurements

Keywords: stereoscopic DIC, elliptic bulge test, ellipsoidal geometry, fibrous connective membrane, liver capsule

Elliptic bulge tests are conducted on liver capsule, a fibrous connective membrane, associated with a field measurement method to assess the global geometry of the samples during the tests. The experimental set up is derived from a previous experimental campaign of bulge tests under microscope. Here, a stereoscopic Digital Image Correlation (DIC) system is used to measure global parameters on the test and investigate some assumptions made on the testing conditions which could not been assessed with microscopic measurements. In particular, the assumption of an ellipsoidal shape of the inflated membrane is tested by comparing the actual sample shape measured by stereoscopic DIC with an idealized ellipsoidal shape.

Results indicate that a rather constant gap exists between the idealized and actual position. The approximation in the calculation of a macroscopic strain through analytical modeling of the test is estimated here. The study of the liver capsule case shows that important differences can be observed in strain calculation depending on the method and assumptions taken. Therefore, analytical modeling of mechanical tests through ellipsoidal approximation needs to be carefully evaluated in every application. Here the field measurement allows assessing the validity of these modeling assumptions. Moreover, it gives precious details about the boundary conditions of the bulge test and revealed the heterogeneous clamping, highlighted by strain concentrations.

Introduction

Experimental protocols coupling microstructural observations to mechanical loading often require taking strong assumptions about the global geometry of the samples during the tests. It comes from the very local nature of the measurements, since it is often very complicated to combine simultaneous macroscopic and microscopic observations. However, in order to use the structural data obtained on a representative volume element (RVE), the local behavior identified has to be homogenized, and to do so, knowledge about the global behavior of the considered material is needed.

Therefore, in the case of bulge testing for instance, assumptions are made about the shape of the deformed sample to compute strain and stress from the measurement of the apex displacement [START_REF] Drexler | Comparison of mechanical behavior among the extrapulmonary arteries from rats[END_REF][START_REF] Yousif | Plastic deformation and failure of thin elliptical diaphragms[END_REF][START_REF] Rees | Plastic flow in the elliptical bulge test[END_REF]. The geometry is thus idealized in order to use the mathematical formulations developed for membranes, shells or plates [START_REF] Flügge | Stresses in Shells[END_REF][START_REF] Feng | Viscoelastic Behavior of Elastomeric Membranes[END_REF].

An alternative to these assumptions is to use a field measurement method such as stereoscopic Digital Image Correlation (DIC). Stereoscopic DIC has often been used in the field of material characterization to perform field measurements at the macroscopic scale, because it allows obtaining information about the whole surface of a sample, and not local measurements of mechanical parameters that are usually obtained through sensors positioned in a specific location. Therefore, with these fields measurements, 2D strain fields can be obtained over the entire surface of the considered material from direct shape measurements and 2D stress fields can be computed numerically without any assumptions on geometrical symmetry [START_REF] Davis | Pointwise characterization of the elastic properties of planar soft tissues: application to ascending thoracic aneurysms[END_REF][START_REF] Genovese | Digital image correlation-based point-wise inverse characterization of heterogeneous material properties of gallbladder in vitro[END_REF]. DIC has thus intensively been used to study various biological materials. For instance, Tonge et al. [START_REF] Tonge | Full-field bulge test for planar anisotropic tissues: part I--experimental methods applied to human skin tissue[END_REF] use a field measurement by stereocorrelation to study the skin anisotropy. Strain fields of the sclera [START_REF] Myers | The inflation response of the posterior bovine sclera[END_REF], liver capsule [START_REF] Brunon | Characterization of the nonlinear behaviour and the failure of human liver capsule through inflation tests[END_REF], parenchyma [START_REF] Brunon | Mechanical characterization of liver capsule through uniaxial quasi-static tensile tests until failure[END_REF][START_REF] Gao | Estimating zero-strain states of very soft tissue under gravity loading using digital image correlation[END_REF] or skin [START_REF] Jacquemoud | Methodology to determine failure characteristics of planar soft tissues using a dynamic tensile test[END_REF] have also been investigated and obtained by DIC.

Therefore, the aim of this study is to use a stereoscopic DIC method to estimate the validity of the assumption of an idealized ellipsoidal deformed shape during an elliptic bulge test conducted on a thin connective membrane. To this end, we chose to adapt the bulge test set up presented in Jayyosi et al. [START_REF] Jayyosi | Characterizing liver capsule microstructure via in situ bulge test coupled with multiphoton imaging[END_REF] to perform a field measurement on the whole sample and thus, getting information on the macroscopic behavior, following the approach developed by Brunon [START_REF] Brunon | Characterization of the nonlinear behaviour and the failure of human liver capsule through inflation tests[END_REF][START_REF] Brunon | Caractérisation et modélisation de la rupture des tissus hépatiques[END_REF]. This global measurement allows checking the hypotheses that are needed to perform an analytical modeling of the tests. These kind of modeling are required when the global field measurement is not accessible, for example in the case of experiments involving a more local imaging method than stereoscopic DIC as reported in Perrini et al. 2013 [16] or Jayyosi et al. 2015 [START_REF] Jayyosi | Characterizing liver capsule microstructure via in situ bulge test coupled with multiphoton imaging[END_REF].

Material and Methods

Samples preparation

The bulge tests have been performed on samples from 2 different human livers (♂ 86 and ♂ 92 years old), which have been harvested at the Department of Anatomy of the University of Rockefeller (DUAR). Informed consent of the subjects and ethical approval have been obtained since livers were harvested from subjects within the French voluntary corpse donation to science program. Following the same method used in Jayyosi et al. [START_REF] Jayyosi | Characterizing liver capsule microstructure via in situ bulge test coupled with multiphoton imaging[END_REF], circular samples of approximately 50mm diameter, have been drilled for the clamping screws in 3 different shape patterns (R1 for circular samples, R2 for elliptic with a ratio 2 between major and minor axis, and R4 for elliptic samples with a ratio of 4). Figure 1a shows the dimensions of the test set up that determined the sizes of the samples. In total, 3 samples have been tested from liver #1 (1xR1 1xR2 and 1xR4) and 8 from liver #2 (3xR1, 3xR2 and 2xR4). Samples from liver #1 were frozen at -20°C then thawed at room temperature before testing.

Samples from liver #2 were tested fresh, i.e. within 3-4 days after the death, while being stored at 4°C on a tissue soaked with PBS to prevent dehydration. DIC measurements are based on the analysis of random patterns on surfaces. Therefore, liver capsules have been spread over a silicone joint to be easily handled (figure 1b), and small drops of paint have been dropped off by gravity by spraying in parallel manner to the plane of the sample (figure 1c). Considering the sizes of the samples, the stains of the random pattern have to be very small to get a pattern heterogeneous enough to allow a good resolution for a precise computation on the whole area of observation. By spraying in a parallel manner to the sample, we can get rid of the large drops of paint that are ejected from the spray, and keep only the small drops of the paint mist, which drop off with gravity (figure 1c). Here, the average sizes of stains are in the order of 0.5-1 mm.

Stereoscopic DIC set up

The experimental set up previously developed for the bulge tests under microscope presented in Jayyosi et al. [START_REF] Jayyosi | Characterizing liver capsule microstructure via in situ bulge test coupled with multiphoton imaging[END_REF] has been reused. Briefly, after preparation, samples on silicone joints are clamped between a top ring and a bottom cylinder thanks to 8 clamping screws equally distributed around the inflation area, and the ensemble is placed into a PMMA (polymethyl methacrylate) inflation base (figure 2). As the liver capsule is translucent, the PMMA inflation base allows positioning the sample between the lighting and the camera. Thus, stereoscopic DIC is performed on the shadow of the random pattern projected on the cameras (figure 2). This particular set up allows avoiding problems of brightness or uniformity of lighting.

Two cameras (Fastcam SA3, Photron, Tokyo, Japan) have been used to record the test, each equipped with a 1024 x 1024 pixel sensor. 105 mm macro objectives have been added to the cameras to adapt the field of view to our samples' dimensions. Depth of field has been set to be able to follow an out of plane displacement of 10 mm, allowing recording the whole test (based on the indicative values obtained from the tests under microscope that indicated an out of plane displacement of 6-8 mm before failure). Cameras have been positioned to perform the correlation on the largest area on the sample, including the edges (figure 2). Then, a calibration has been performed thanks to a dedicated target to compute the relative position of the cameras, needed for the stereoscopic calculation.

The image acquisition has been performed at 50 Hz until failure to assess the rupture creation/propagation. Pressure was acquired at 100 Hz with accuracy of ±1 mbar. As it was the case for the bulge tests under microscope, inflation is performed with compressed air. Loading is conducted with regular pressure increments up to failure, thanks to a manual valve.

Data processing

Displacement and Green-Lagrange strain fields are computed by the image correlation software VIC3D (Correlated Solutions, Inc.) [START_REF] Sutton | Image correlation for shape, motion and deformation measurements[END_REF]. Correlation is performed on a subset of 500-600 images equally distributed among all the images taken during the test, to keep a reasonable computation time. On these image sequences, the displacements are calculated on subsets of 31 x 31 pixel wide and spaced of 7 pixels from one to another (figure 3a). A rotation of the computation basis has been applied to match with the basis linked to the sample (X axis on the major axis of the ellipse, Y axis on the minor axis of the ellipse and the Z axis corresponds to the vertical).

An analytical modeling of the inflated geometry is made for each test based on an approximation by an ellipsoid (equation 1). The dimensions of the approximated ellipsoid correspond to the dimensions of the ideal boundary conditions in the X Y plane (R1, R2 or R4), associated with the vertical position of the apex.

with θ ϵ [0 , π] and φ ϵ [-π , π], and where a and b are the dimension of the clamping in the XY plane (figure 1a), and c the vertical position of the apex.

In order to test the validity of this approximation of the inflated dome geometry, the error between the real and approximated position of the capsule has been calculated, performing an orthogonal projection of the DIC measured surface on the idealized ellipsoidal surface. To do so, for each i point of the DIC surface (~2000-3000 points depending on the sample type), we look for the minimum distance di with every j points of the analytical surface (40000 points).

Then a mean error is calculated at each pressure level as an average of the distance errors of each surface point.

In order to compare it to the DIC strain, macroscopic strain derived from the deformed ellipsoid (referred to as ellipsoidal strain in the following) is also computed based on that ellipsoid geometry assumption. The ellipsoid meridian length in X and Y directions are calculated from the second formula of Ramanujan [START_REF] Ramanujan | Modular equations and approximations to pi[END_REF], which gives for the X direction:

where a and c are the previously mentioned dimensions of the ellipsoid (a is replaced by b in this formula to calculate the meridian length in the Y direction).

Then, the stretch ratios in X and Y directions are given by: And Green-Lagrange strain at the apex in X and Y directions is calculated as follow:

� 𝑥𝑥 = 𝑎𝑎 cos 𝜃𝜃 cos 𝜑𝜑 𝑦𝑦 = 𝑏𝑏 cos 𝜃𝜃 sin 𝜑𝜑 𝑧𝑧 = 𝑐𝑐 sin 𝜃𝜃 (1)
𝑑𝑑 𝑖𝑖 = min 𝑗𝑗 � �𝑋𝑋 𝐷𝐷𝐷𝐷𝐷𝐷 𝑖𝑖 -𝑋𝑋 𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴 𝑗𝑗 � 2 + �𝑌𝑌 𝐷𝐷𝐷𝐷𝐷𝐷 𝑖𝑖 -𝑌𝑌 𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴 𝑗𝑗 � 2 + �𝑍𝑍 𝐷𝐷𝐷𝐷𝐷𝐷 𝑖𝑖 -𝑍𝑍 𝐴𝐴𝐴𝐴𝐴𝐴𝐴𝐴 𝑗𝑗 � 2 (2)
𝐿𝐿 𝑥𝑥 ≈ 𝜋𝜋(𝑐𝑐 + 𝑎𝑎) 4 �1 + 3ℎ 10 + √4 -3ℎ � (3)
where

ℎ = � 𝑐𝑐 -𝑎𝑎 𝑐𝑐 + 𝑎𝑎 � 2 (4)
� 𝜆𝜆 𝑥𝑥 = 𝐿𝐿 𝑥𝑥 𝑎𝑎 𝜆𝜆 𝑦𝑦 = 𝐿𝐿 𝑦𝑦 𝑏𝑏 (5)
� 𝐸𝐸 𝑥𝑥𝑥𝑥 = 1 2 �𝜆𝜆 𝑥𝑥 2 -1� 𝐸𝐸 𝑦𝑦𝑦𝑦 = 1 2 �𝜆𝜆 𝑦𝑦 2 -1� (6)

Results

Evolution of the apex vertical displacement with pressure increase

The initial state for image correlation has been chosen to correspond with the state where the capsule presents a uniform curvature, which is achieved for a pressure of 0.02 bar, in accordance with what was previously reported by Brunon [15]. This initial pressure of 0.02 bar is almost 20 times smaller compared to rupture pressure.

An example of the global shape of the inflated capsule is given on figure 4. We notice that the correlation does not work systematically on the edges near the clamping due to small droplets of saline solution which can gather up in these areas and erase the random pattern locally, making the correlation fail. However, the particular positioning of the cameras allows keeping a large area for the correlation computation, as the left and right cameras show a large common area in their respective field of view (figure 4).

The evolution of the vertical displacement of the apex with the pressure increase presents some variability depending on the sample considered. Figure 5 shows the apex vertical displacementpressure curves. The sudden increase in vertical displacement for sample #1 at the end of the test comes from slippage of the capsule in the clamping. Globally speaking, the specific geometry of the different test types induces a greater vertical displacement for the circular samples than for the elliptic ones.

Geometry approximation by an ellipsoid

Figure 6a illustrates the evolution of the mean error between the shape measured by stereoscopic DIC, and the idealized ellipsoidal shape. This error is within the range of 0.25-0.4 mm and does not vary much with the level of pressure (except for a few cases). Figure 6b gives an example of the standard deviation associated to the mean errors computed over the all surface area for a sample of each type R1, R2 and R4. The standard deviation appears quite important as the error between the position of the two surfaces varies depending on the location on the ellipsoid (it is null at the apex and increase as we get closer from the outer perimeter). An example of the comparison between idealized and real shapes is given on figure 7.

Stereoscopic DIC strain versus ellipsoidal strain

Post processing reveals a mean deformation rate around 2.10 -3 s -1 in the Y direction (which correspond to the direction of the minor axis of the ellipses).

Figure 8 presents the comparison between the stereoscopic DIC strain and the ellipsoidal strain at the apex in the X and Y directions, for all the samples tested. A distinction is made between the different sample types R1, R2 and R4. An important gap can be observed between the two calculated strains over the samples considered, particularly for the R4 samples. For the R1 samples the stereo and ellipsoidal strains match properly except for the first sample (which experienced sliding at the end as seen on figure 5).

Strain fields

Figure 9 shows an example of the Green-Lagrange strain fields Exx and Eyy as well as principal strain E1 and E2 for each sample type. Strain fields often show some strain concentrations around the boundaries at the clamping. The non-uniformity of the tightening of the clamping screws induces heterogeneities in the strain fields. However, these variations are successfully assessed by the field measurements by stereoscopic DIC, which allows getting a realistic idea of the effective strain state of the sample. It therefore gives an accurate assessment of the real boundary conditions of the bulge test. We notice that the directions of principal strain for the elliptic samples match properly the X and Y directions of the ellipse when the tightening is uniform (figure 9b). On the other hand, these directions can differ from the directions of the ellipse when the tightening induces heterogeneities in the strain fields (figure 9c).

Discussion

This experimental campaign aim was to quantify at the macroscopic scale the loading conditions applied during the experiments under microscope described in Jayyosi et al. [START_REF] Jayyosi | Characterizing liver capsule microstructure via in situ bulge test coupled with multiphoton imaging[END_REF]. The displacement field measurements on the whole sample allowed assessing the real sample geometry and boundary conditions during the bulge test.

The relative error between the real position of the capsule and the position approximated by an ellipsoid indicates that the analytical modeling of the test must be carried out carefully. Using ellipsoid parametric equations, the stretch ratios can be computed to estimate a global strain measure on the different tests from the apex displacement value solely. However, the error on the position of the capsule made with the ellipsoidal assumption induces uncertainty on the macroscopic strain calculation based on such analytical modeling. Indeed, a mean error of 0.35 mm on the apex position, which is globally observed on the samples tested here as seen on figure 6, should create an uncertainty between 0-5 % on the macroscopic strain value calculated. Therefore, the estimation of macroscopic strain from that perspective remains an approximation.

In fact, the difference between stereoscopic strain and ellipsoidal strain exceed those 5 % strain to reach 10-15% in extreme cases as seen on figure 8. The additional assumptions taken to calculate strain (especially in the stretch ratios calculation) combined with the uncertainty on the apex position lead to significant errors on ellipsoidal strain estimation.

It is noteworthy that here, the ellipsoidal approximation made is not optimal as the DIC surface has not been fit with an ellipsoid with the objective of minimizing the error between the two surfaces. Figure 6b highlights this point by showing the important standard deviation of errors over the all surface. The aim was to compare the measured DIC shape with an ellipsoidal approximation made solely from the knowledge of the apex displacement and the boundary conditions, to evaluate the analytical modeling of the test based on these scarce measurements. A more accurate analytical modeling of the test could be performed from the data of the field measurements by looking for an optimal ellipsoidal fit as in [START_REF] Tonge | Full-field bulge test for planar anisotropic tissues: part I--experimental methods applied to human skin tissue[END_REF][START_REF] Brunon | Characterization of the nonlinear behaviour and the failure of human liver capsule through inflation tests[END_REF].

Another source of explanation for the difference observed between the deformed shape and the analytical shape could come from the heterogeneity or the anisotropy of the material. Indeed, Tonge et al. [START_REF] Tonge | Full-field bulge test for planar anisotropic tissues: part I--experimental methods applied to human skin tissue[END_REF] showed that skin anisotropy induce an ellipsoidal shape for the inflated sample with circular boundary conditions instead of a spherical shape. Such effect cannot be predicted with the kind of analytical modeling presented here. However, for the case of liver capsule, this effect is not likely to have a significant impact as the microstructural organization of this tissue does not present a clear marked anisotropy [START_REF] Jayyosi | Characterizing liver capsule microstructure via in situ bulge test coupled with multiphoton imaging[END_REF], and the optimized spherical fit of a bulge test of liver capsule from DIC measurements gave satisfying results [START_REF] Brunon | Characterization of the nonlinear behaviour and the failure of human liver capsule through inflation tests[END_REF].

A high variability of behavior among the different samples is observed, as seen on figure 5. It highlights the difficulty of identifying mechanical properties for biological tissues. Indeed, the intersamples variability within one liver is already important. To that variability, it must be added the dispersion that comes from the inter-individual variability, which makes it even more difficult to identify a unique law to describe the mechanical behavior. A part of the variability observed here also comes from the specific testing conditions. The loading rate is not completely identical from one sample to another, despite the experimental care to make every test repeatable. Thus, the tests comparison is made without any considerations to viscous effects.

Moreover, the clamping realized by screws equally distributed around the sample is not strictly uniform, and probably has an effect on the displacements measured by stereoscopic DIC (this is the case for sample #1 for example that slipped at the end of the test). This non uniform clamping induces a non-homogeneous deformation of the silicone joint which creates the heterogeneities in the boundary conditions. The problem mentioned in the Brunon thesis [START_REF] Brunon | Caractérisation et modélisation de la rupture des tissus hépatiques[END_REF], about the boundary conditions control by realizing a precise and uniform clamping has not entirely been solved here. Positioning more clamping screws allowed limiting slippage of the sample, by handling more precisely the boundary conditions. However, the field measurement still reveals the existence of heterogeneous strain concentrations on the borders of the samples, as highlighted by figure 9. These stress concentrations in the clamping screws question one of the perks of bulge testing, which is to load the sample far from the edges, to avoid rupture initiation due to cutting. Still, it is noteworthy that the field measurement allows assessing and quantifying this phenomenon, and by doing so, knowing the actual strain state of the material during the test.

The heterogeneity of the strain fields also points out a limitation of the ellipsoidal analytical modeling of the test since these disparities in strain cannot be accounted for with such modeling. Therefore, the gap between stereoscopic DIC strain and ellipsoidal strain, seen on figure 8, comes also from this phenomenon. This kind of analytical modeling based on scarce measurements does not take into account the heterogeneous boundary conditions. In that regard, field measurements appear as a valuable and practical solution when a detailed knowledge of the all sample mechanical state is needed. However, when such methods as stereoscopic DIC cannot be used, the alternative methods and hypotheses taken to infer information from scarcer measurements need be to carefully investigated to get a clear idea on how they impact the obtained results [START_REF] Murienne | A comparison of 2D and 3D digital image correlation for a membrane under inflation[END_REF].

Conclusion

The implementation of a macroscopic field measurements coupled to the bulged tests, previously developed for microscopic tests, allowed accurately assessing the experimental conditions of the inflation test. Especially, the global geometry of the inflated dome and the boundary conditions during the test have been measured, thanks to the stereoscopic DIC method.

Stress concentrations occurred on the edges during the inflation, thus revealing a limitation of the experimental protocol, associated with the non-uniform clamping of the sample.

The field measurement by stereoscopic DIC offers the possibility to evaluate the performance of the ellipsoidal shape assumption. This assumption is necessary to compute macroscopic strain from scarce measurements, which is the case in the microscopic tests performed under microscope, where just the measure of the apex displacement is made. The evaluation of this analytical modeling performance provides knowledge about the accuracy of this macroscopic strain estimation. Eventually, the comparison between macroscopic and microscopic strains will allow linking the microstructural properties to the macroscopic behavior of the liver capsule.

Figure 1 :

 1 Figure 1: a) Dimensions of the bulge test set up: R1: a=b=10mm, R2: a=12mm and b=6mm, R4: a=20mm and b=5mm. b) Samples spread over the silicone joint on which we perform the random pattern for stereoscopic DIC. c) Random pattern realization technique with drop off of paint drops by gravity.

Figure 2 :

 2 Figure 2: Lay out of the experimental set up with the PMMA inflation base positioned between the lighting spots and the 2 cameras for stereoscopic correlation.

Figure 3 :Figure 4 :

 34 Figure 3: a) Initial subset of the image showing the random pattern before deformation. b) Deformed subset with updated position (VIC3D®)

Figure 5 :

 5 Figure 5: Evolution of vertical displacement of the liver capsule apex with pressure increase, for the different types of loading considered (R1, R2 and R4).

Figure 6 :

 6 Figure 6: a) Mean error between the actual position measured by stereoscopic DIC and the approximated position of the capsule by an ellipsoid. The approximated ellipsoid is based on the theoretical dimensions of the clamping and the vertical displacement of the apex. b) Standard deviation associated to the error calculated over the all surfaces for sample #3 (R1), #7(R2) and #9 (R4).

Figure 7 :

 7 Figure 7: Geometry of the inflated capsule for a R4 sample (blue) put in relation with the ellipsoidal approximation (red).

Figure 8 :

 8 Figure 8:Comparison between stereoscopic DIC strain (-) and ellipsoidal strain (..) in the X direction (left column) and Y direction (right column), for all the samples tested (sorted out between R1, R2 and R4).

Figure 9 :

 9 Figure 9: Example of Green-Lagrange strain fields in the X, Y and principal directions for a R1 sample (a), R2 sample (b) and R4 sample (c).

Fig

  Fig 1.a

Fig 2 13

 13 Fig 2

Fig. 3

 3 Fig. 3.a

Fig. 6

 6 Fig. 6.a

Fig. 8

 8 Fig.8

Fig. 9

 9 Fig. 9.a

Acknowledgments

This work was supported by the Programme Avenir Lyon Saint-Etienne (ANR-11-IDEX-0007) of Université de Lyon, within the program "Investissements d'Avenir" operated by the French National Research Agency (ANR).

Disclosure

None of the authors have any professional or financial conflict of interest.


Eugene Asarin

Thierry Cachat

Alexander Seliverstov

Tayssir Touili
email: touili@liafa.jussieu.fr

Vassily Lyubetsky
email: lyubetsk@iitp.ru

Attenuation Regulation as a Term Rewriting System ⋆

The classical attenuation regulation of gene expression in bacteria is considered. We propose to represent the secondary RNA structure in the leader region of a gene or an operon by a term, and we give a probabilistic term rewriting system modeling the whole process of such a regulation. ⋆ The support of CNRS-RAS cooperation agreement 19122 Evolver is gratefully acknowledged. * → means several similar rewritings; repeated window positions (e.g. repetitions of 40, 51 )are replaced by a • symbol; ⊥ means termination. There are 24 helices, denoted by letters from a to x.

Introduction

Modeling the mechanisms of regulation of gene expression, allowing prediction of quantitative characteristics of this expression (such as estimation of the level of expression and concentration of the substrate) is an important research challenge. In a previous work [LRSP06,LPRS07], a model of one particular kind of regulation, the classical attenuation regulation, has been suggested. In that model, the evolution of the secondary RNA structure in the leader region of a gene, and the progress of the ribosome and the polymerase along the RNA/DNA strands, are represented by a very special, elaborated in detail, Markov chain. In this chain the transition probability corresponding to the progress of the ribosome depends on a "control variable" -the concentration of charged tRNA molecules in the cell. All the other probabilities do not depend on the control variable, they can be determined from energy-based considerations. Termination and antitermination (of gene expression) correspond to particular random events in the Markov chain. In [LRSP06], a Monte-Carlo simulation of this Markov chain led to biologically realistic dependence of termination probability from the control variable. Due to a large size and a complex structure of the Markov chain, its simulation is a heavy computational task, but it was successfully solved, and a software tool called Rnamodel simulates one trajectory in fractions of a second [LRSP06,RNA]. However, the approach based on the direct description of the Markov chain and its simulation has some limitations, especially for a theoretical analysis. Biologically, it would be nice to have a more structured and compact representation of the Markov chain and its instantaneous probability distributions over all states at every instant, or only for sufficiently large time, or only probabilities of the two biologically important events -termination and antitermination.

Note that the problem of modeling the classical attenuation regulation, as stated in [LRSP06] and in the current article, is related to the representation of the transient behavior of the secondary structure on a sliding window on the RNA strand between the ribosome and the polymerase (see below for details). This differs from the kinetics of the secondary RNA structure on a fixed nucleotide sequence for unlimited time, i.e. unlimited number of steps, investigated in many papers. The structure that appears after a large amount of time is called equilibrium secondary RNA structure, it corresponds to a minimum of energy, see e.g. [Zuk03,FFHS00]. The tool Rnamodel has also the function of determining this equilibrium structure and its energy as a special part of the full model in [LRSP06]. However, real structures that appear on the RNA strand during the regulation process are far from the equilibrium and their energies are far from minimal.

In this article we discover a regular internal structure of the Markov chain describing the classical attenuation regulation. We show that it can be represented as a probabilistic term rewriting system for a particular type of terms. The set of rewriting rules can be large, but all of them are generated by a small set of (five) metarules. In fact we give the full description of the metarules and explain how to generate all the rules for the case of classical attenuation regulation.

Potential benefits of such a representation are multiple:

-easier and more precise modeling of regulation mechanisms depending on the dynamics of the secondary structure; -compact description of such mechanisms, perhaps in dedicated languages, and hence a better biological understanding of regulation processes; -convenient representation of secondary structures by terms; -specific analysis and simulation methods for rewriting systems. This article is structured as follows. In section 2 we describe shortly the biological phenomenon that we want to model: the mechanism of classical attenuation regulation (CAR). In section 3 we introduce a class of terms and probabilistic term rewriting systems. In section 4 we represent a qualitative metamodel of the biological mechanism of CAR by a term rewriting system. In section 5 we refine the previous system and decorate its transitions with rates, thus obtaining a representation of the Markov chain by a probabilistic term rewriting system. In section 6 we show some simulation results. In section 7 we discuss some related work on term rewriting and its applications. In section 8 we conclude with a discussion of perspectives of the rewriting approach to modeling the mechanisms involving RNA secondary structures, especially regulation.

Classical attenuation regulation

To begin with, we recall some well-known biological facts about the biological phenomenon playing the central role in this article.

The expression of a group of structural genes (that is synthesis of the corresponding proteins, which are ferments for a chemical reaction) can be regulated by a sequence of nucleotides placed on the DNA upstream inside the so called leader region of the genes [SB91]. This subsequence of the leader region is called the regulatory region. In this article we deal with one particular type of regulation, classical attenuation regulation (CAR) in bacteria. This regulation mechanism concerns structural genes (groups of genes -operons) that produce proteins which catalyze the synthesis of amino acids. The classical attenuation allows to activate such an operon when the cell contains a small concentration of the amino acid, to deactivate the operon whenever this concentration increases, and to do it fast. The mechanism of CAR involves several actors: the regulatory region on the DNA, its copy on the RNA, the ribosome, and a ferment called RNA polymerase (see Fig. 1).

Rib Pol

str. genes

DNA RNA Q Q Q' Q' Q" Q"
Figure1. Classical attenuation regulation. The RNA polymerase Pol transcribes the regulatory region Q, the ribosome Rib translates the leader peptide gene Q ′ . The movement of Rib on regulatory codons Q ′′ is controlled by the concentration of charged tRNA. The secondary RNA structure ω between Rib and Pol brakes Pol and pushes it off the chain. If Pol reaches the structural genes, then they are expressed, i.e. transcribed and then translated. Note that in both the DNA and the RNA, we use Q, Q ′ and Q ′′ to denote the regulatory region, the leader peptide gene, and the regulatory codons, respectively.

For structural genes to be expressed two concurrent processes should succeed: the regulatory region Q should be transcribed creating an RNA by RNA polymerase. At the same time the ribosome should be bound to the very beginning of the freshly created segment Q ′ (called the leader peptide gene) in the regulatory region Q on the RNA and starts translation of this leader peptide gene to an auxiliary protein. The essential part of the regulation process takes place when the ribosome moves on Q ′ on the RNA and the polymerase moves somewhere downstream of the ribosome on Q on the DNA.

The ribosome moves "rightwards" (formally speaking, in the direction from the 5 ′ to the 3 ′ end) on a segment Q ′ of the sequence Q. Its speed is constant except on a subsequence Q ′′ (regulatory codons) where it depends directly on the concentration of the amino acid (via charged tRNA concentration). To the right of the ribosome and independently of it, the polymerase moves rightwards on Q. Between the ribosome and the polymerase a secondary structure ω is formed on the RNA. This structure consists in pairing of some nucleotides, and it changes very fast. An important effect of the secondary structure ω consists in slowing down the movement of the polymerase. There are two possible scenarios:

-When ω is strong enough, its "braking" action on the polymerase increases, and moreover, the polymerase can slip off the DNA (this can only happen on so-called T-rich sequence, where the connection of the polymerase and the DNA weakens). Such an event is called termination, and in this case the structural genes are not expressed: the transcription of the regulatory region is aborted, the structural genes are not transcribed and therefore not translated.

-Another possibility is that the ribosome moves fast enough to weaken or partly destroy most of the structure ω. In this case the polymerase safely traverses the T-rich sequence, and arrives to the end of the leader region Q.

Next, the polymerase enters the structural genes, and their transcription, followed by translation are unavoidable. This event is called antitermination and in this case the structural genes are expressed.

In the rest of this article we build a qualitative and a quantitative models of the regulation process described above.

3 Terms and rewriting systems

Unranked unordered terms

Let Σ be a finite set of function symbols and X an enumerable set of variables (standing for sets of terms). The set T Σ [X ] of terms over Σ and X is the smallest set that satisfies:

-Σ ⊆ T Σ [X ], -{f (x) | f ∈ Σ ∧ x ∈ X } ⊆ T Σ [X ], -if f ∈ Σ and s ⊆ T Σ [X ] is a set of terms, then f (s) is in T Σ [X ].
By definition we also put f (∅) = f for f ∈ Σ. For convenience we write f (g, h(e)) instead of f ({g, h({e})}). However one should remember that the coma-separated terms are unordered.

Example 1. Let Σ = {e, f, g, h} and X = {x, y, z, . . . }, then the followings are terms in T Σ [X ]: f (g, h(e)), f (f (x)) and e(g, f ).

Note that we consider function symbols of variable arity. T Σ stands for T Σ [∅]. Terms in T Σ are called ground terms. Variables are used only to define substitution and rewriting rules. The "real" terms are ground terms. A substitution σ is a mapping from X to 2 TΣ [X ] , written as σ = {x 1 → T 1 , . . . , x n → T n }, where T i , 1 ≤ i ≤ n, is a finite set of terms that substitutes the variable x i . The term obtained by applying the substitution σ to a term t is written tσ. We call it an instance of t.

Let R be a rule of the form l → r, where l and r are terms in T Σ [X ]. For ground terms t, t ′ we write t → R t ′ if there exists a substitution σ such that t ′ can be obtained from t by replacing an occurrence of the subterm lσ by rσ. → R defines a relation between ground terms. Let → * R be the reflexive transitive closure of → R .

Example 2. Let R = l → r with l = f (x, e), r = f (g(x), e) and t = e(f (h, e)), then t → R t ′ where t ′ = e(f (g(h), e)).

A term rewriting system (TRS) is a finite set of rules of the form l → r. Given a TRS R and a set of terms I ⊂ T Σ , the language R * (I) is defined as the set of all ground terms that can be obtained from the terms in I by applying a finite number of times the rules from R, i.e., R *

(I) = {t ∈ T Σ | ∃t ′ ∈ I, t ′ → * R t}. Example 3. Let R = {f (x) → g(f (x))} and I = {f (e, h)}, then R * (I) = {g n (f (e, h)) | n ∈ IN}.

Probabilistic Term Rewriting Systems

A Continuous Time Markov Chain is a pair (S, ρ), where S is a finite or enumerable set of states and ρ : S × S → [0, ∞) is the rate matrix. For s, s ′ ∈ S, ρ(s, s ′ ) > 0 means that there is a transition between states s and s ′ , and that the probability for moving from s to s ′ within t time units is equal to 1e -ρ(s,s ′ )•t . If a state s has more than one outgoing transition (i.e., if there exist more than one state s ′ for which ρ(s, s ′ ) > 0) there exists a race between these transitions and the probability for moving from s to s ′ within t time units is equal to

ρ(s,s ′ ) E(s) 1 -e -E(s)•t , where E(s) = s ′ ∈S ρ(s, s ′ ).
A (continuous time) Probabilistic term rewriting system (PTRS) over Σ ∪ X is a (finite) set of rules of the form l Λ --→ r, where l and r are terms in T Σ [X ], and Λ ∈ (0, ∞) is a rate.

A PTRS R over Σ ∪ X defines a continuous time Markov chain on ground terms M = (T Σ , ρ), where ρ(t, t ′ ) = Λ iff there exists a rule l Λ --→ r ∈ R such that t → R t ′ , where R is the "non probabilistic" rule l → r.

Remark 1. If there are several rules (or several instances of the same rule) that lead from t to t ′ , then ρ(t, t ′ ) = Λ, where the sum is taken over all such rules or instances.

Metamodel

We want to model the phenomenon of the classical attenuation regulation described in section 2.

We suppose that a regulatory region Q (see Fig. 1) is given and fixed in the sequel, it is a sequence (word) Q ∈ {A, C, G, T} * , the letters of this alphabet are called nucleotides. We denote by |x| the length of any word x and x i the ith letter of x, so x = x 1 x 2 . . . x |x| . The sequence Q can be folded1 in a way that some nucleotides of Q are paired: A with T and C with G. The complement of a nucleotide is written using a bar: A = T, T = A, C = G, G = C. We look in Q for subwords ("stems") of the form

Q A Q A+1 . . . Q B and Q C Q C+1 . . . Q D such that B -A = D -C, A + 3 ≤ B, B + 3 ≤ C
(1)

Q A = Q D , Q A+1 = Q D-1 , . . . Q B = Q C .
Any pair of such stems forms a hypohelix (see Figure 2, where the labels A i , B i , C i and D i are positions in the word Q).

1 A1 B1 C1 D1 |Q| f Figure2. One hypohelix f .
We describe a hypohelix f by a tuple of its stems' extremities f = (A, B, C, D), and we introduce the following notations:

stem(f ) = [A, B] ∪ [C, D], loop(f ) = [B + 1, C -1], supp(f ) = [A, D].
There is a ribosome at some position on Q ′ and an RNA polymerase somewhere to the right of it. Both move to the right, in one step the ribosome moves by three successive nucleotides and the polymerase by one nucleotide. The window w = (R, P ) represents the segment of RNA from the first position R after the end of the ribosome to the last position P before the beginning of the polymerase. In fact the folding of the RNA sequence Q can only happen within the current window, i.e. between positions R and P . When the ribosome advances to the right, it can destroy the leftmost hypohelix of a current configuration, because it consumes the first three letters of the window. On the other hand any polymerase move adds one new letter to the window.

Formally a window has the form w = (R, P ) with R, P ∈ IN. The following constraints should be satisfied:

13 ≤ R ≤ P ≤ |Q| (2)
Thus, the window is moving and changing its length.

Let W = {w = (R, P ) | conditions (2) are satisfied} be the alphabet of all windows. We define

stem(w) = ∅, loop(w) = [R, P ], supp(w) = [R, P ] .
We will write terms over the alphabet Σ of all hypohelices and all windows:

Σ = H ∪ W where H = {f = (A, B, C, D) | conditions (1) are satisfied}.
We consider only terms of the form w(. . . ) for some w ∈ W (rooted by some window w). According to the conditions that we will define next, a symbol f = (A, B, C, D) can appear in a term w(. . . ) only if R ≤ A and D ≤ P , where w = (R, P ).

We say that a hypohelix f is embedded in g (which can be a hypohelix or a window), written f ≺ g, if supp(f ) ⊆ loop(g). Two hypohelices f and g are disjoint, written f ⊲⊳ g, if supp(f ) ∩ supp(g) = ∅. We call f and g unknotted if either one of them is embedded in the other or they are disjoint. We say that g

= (A 2 , B 2 , C 2 , D 2 ) is an extension of f = (A 1 , B 1 , C 1 , D 1 ), denoted f ⊑ g, if [A 1 , B 1 ] ⊆ [A 2 , B 2 ] and B 2 -B 1 = C 1 -C 2 , hence [C 1 , D 1 ] ⊆ [C 2 , D 2 ],
and the pairing in g is an extension of that in f . See Figure 3.

We call a term t over Σ well-formed if it satisfies the following conditions:

(compatibility) any f and g appearing in t are unknotted, in particular any f can appear at most once, (ordering) if f and g occur in t, then f ≺ g iff f is in the scope of g.

The combination of two hypohelices in Figure 4 is biologically feasible, but according to our rules these hypohelices are incompatible. We believe that this restriction (crucial for representation by terms) does not undermine significantly the accuracy of the model.

Notice, that a well-formed term of the form w(. . . ) (rooted by some window w) contains only hypohelices from

Σ w = {f ∈ H | f ≺ w}.
This simple observation greatly simplifies the simulation process.

In [LRSP06] an additional maximality condition is imposed. Using the terminology of this article, it requires that no hypohelix f in t can be replaced by its proper extension without creating an overlapping. Here we do not impose this restriction.

Each well-formed term represents a possible secondary RNA structure in a window in Q: the set of hypohelices that are present in this window. It could be

A1 B1 C1 D1 f A2 B2 C2 D2 g A2 B2 C2 D2 f A2 B2 C2 D2 g
Figure3. Relative positions of two hypohelices f and g: f ≺ g and f ⊲⊳ g.

Here f = (A 1 , B 1 , C 1 , D 1 ) and g = (A 2 , B 2 , C 2 , D 2 ). On the left B 2 < A 1 and D 1 < C 2 , on the right D 2 < A 2 . A1 B1 C1 D1 A2 B2 C2 D2 Figure4. Pseudo-knot: A 1 < B 1 < A 2 < B 2 < C 1 < D 1 < C 2 < D 2 .
Such configurations are not allowed in our model.

possible to allow knotted hypohelices, and hypohelices of length less than 3, but here we do not consider them. We extend the definitions of ⊲⊳ and ≺: let f be a term and c a set of terms,

c ⊲⊳ f iff ∀g ∈ c (g ⊲⊳ f ), c ≺ f iff ∀g ∈ c (g ≺ f ).
In the former case we say that f and c are disjoint, in the latter that c is embedded into f . We start from a sequence Q without any pairing of nucleotides, this structure is described by a term w() -"an empty window", where w = (13, 13). Our aim is to represent the evolution of the secondary structure in the window, as well as the progress of the ribosome and the polymerase, through rewriting terms starting from w(). Our rewriting system will generate only well-formed terms.

On the whole, there are five rewriting Meta-rules:

-Binding and decomposition of a hypohelix f :

ω = g(c, d) ←→ ω ′ = g(c, f (d)) with c ⊲⊳ f, d ≺ f, f ≺ g, ( 3
)
where c and d are sequences of terms. The concrete rewriting rules -and their rates -depend on c and d, as explained below. -Extension and reduction of a hypohelix

ω = f ←→ ω ′ = g with f ⊑ g. (4)
-The window movement can be described by the following rules, where w = (R, P ):

(R, P )(ω) -→ (R + 3, P )(ω ′ ) , (5) (R, P )(ω) -→ (R, P + 1)(ω) , ( 6)

w(ω) -→ ⊥ . ( 7
)
In the last rule, ⊥ is a special symbol denoting termination. Rules (5) describe the movement of the ribosome. In these rules, ω ′ is obtained from ω by removing only the possible symbol that is incompatible with the new window (R + 3, P ), or replacing it by a "shorter" hypohelix. Indeed, if the leftmost hypohelix in ω starts at a position between R and R + 3, then the movement of the ribosome by three positions to the right will destroy this hypohelix. More formally, if ω ≺ (R + 3, P ), then ω ′ = ω. Otherwise the ribosome destroys the leftmost hypohelix. In this case, there is a single symbol f in ω such that f ≺ (R + 3, P ). Suppose the subterm rooted by f is f (c). Then, ω ′ is obtained by replacing in ω f (c) by either f ′ (c) or c, depending on the size of f , where f ′ ⊑ f . Rules 6 describe the movement of the polymerase. Note that if the polymerase reaches a position P + 1 where the structural genes are expressed, then we reach antitermination and the gene is expressed. Now, we introduce the rates of the five rewriting rules.

Let h(f 1 ( * ), . . . , f n ( * )) be a term. Then the free loop length of the hypohelix h in this term is

l h = | loop(h)| - n i=1 | supp(f i )| .
This numeric characteristic corresponds to the number of nucleotides in the loop of the hypohelix h that do not participate in inner hypohelices.

In order to define the rate, we have to consider the concrete rule corresponding to the Metarule (3). For any f, g, c = c 1 (x 1 ), . . . , c m (x m ) and

d = d 1 (y 1 ), . . . , d n (y n ) such that c ⊲⊳ f, d ≺ f, f ≺ g there is a concrete rule ω = g(c 1 (x 1 ), . . . , c m (x m ), d 1 (y 1 ), . . . , d n (y n )) ←→ ω ′ = g(c 1 (x 1 ), . . . , c m (x m ), f (d 1 (y 1 ), . . . , d n (y n ))) (8)
Recall that the subterms are unordered. Similarly the concrete rule corresponding to ( 4) is

ω = a(c 1 (x 1 ), . . . , c m (x m ), f (d 1 (y 1 ), . . . , d n (y n ))) ←→ ω ′ = a(c 1 (x 1 ), . . . , c m (x m ), g(d 1 (y 1 ), . . . , d n (y n ))) (9)
Note that this transformation can change the free loop length of the hypohelix a. The rate of the rules (8-9) is denoted K(ω → ω ′ ), given by

K(ω → ω ′ ) = κ • exp 1 2 (E(ω) -E(ω ′ )) , (10)
where the energy E(ω) = G hel (ω)+G loop (ω), κ is a parameter -usually κ = 10 3 -and

G hel (ω) = 1 RT • h E h and G loop (ω) = h 1.77 • ln(l h + 1) + B , (11)
and h varies over all hypohelices from ω. E h represents the total stacking energy along the hypohelix h. It is the sum of stacking bond energies of the adjacent base pairs of h. B can take three different values depending on the three possible types of the loop of the hypohelix g: terminal loop, single-strand bulge and doublestrand bulge.

A codon is a triple of successive nucleotides. For a sequence Q ′ , each codon is fixed to be either regulatory or non-regulatory. Analogously, each nucleotide in Q is fixed to be either non T-rich or T-rich [LRSP06]. Let s 0 be the "radius" of a ribosome -distance from P-site to the end of the ribosome -usually s 0 = 12, and let s 1 be the "radius" of a polymerase -distance from the 5 ′ end of a polymerase to its transcription center -usually s 1 = 9. The rate of the rule (5) is denoted λ rib and is constant when Rs 0 is a position of a nonregulatory codon, and otherwise λ rib depends on an external parameter c -the concentration of charged tRNA [SB91]. The rate of the rule ( 6) is denoted ν and depends on secondary structure ω in the window. The rule (7) applies only when P + s 1 is a position of a T-rich nucleotide and its rate is denoted µ.

In [LRSP06] the rate of the rule (5) was denoted λ rib and

λ rib (c) = 45 c 1 + c . ( 12
)
The rate of the rule (6) was denoted ν and

ν = 40 -F (ω) . ( 13
)
The rate of the rule (7) was denoted µ and

µ = 1 4 F (ω) . (14)
The function F (ω) in (13-14) for ω = f 1 ( * ), . . . , f n ( * ) depends only on functional symbols (hypohelices) f 1 , . . . , f n , and not on the structure of their arguments denoted by * . More precisely F (ω) = max i F (f i ), where

F (f ) = δ • exp -r(f ) r0 (L 2 ) 2 • (p(f ) -p 0 ) 2 + 1 , ( 15
) with p(f ) ≈ π | supp(f )|
, and r(f ) the "free distance" from f to the end P of the window: for f = (A, B, C, D) and w = (R, P ), we have

r(f ) = R -D - i | supp(f i )| . (16)
Other symbols in equation ( 15) denote constants:

r 0 = 1, δ = 30, L 2 = 27.1, p 0 = 0.18, see [LRSP06].
Note that the rates of the rules depend only on the local configuration as explained above and not on the outside context. In particular it does not depend on instantiations of x 1 , . . . , x m , y 1 , . . . , y n .

Simulation results

atgaaagcaattttcgtactgaaaggttggtggcgcacttcctgaaacgggcagtgt attcaccatgcgtaaagcaatcagatacccagcccgcctaatgagcgggcttttttttg

Figure5. A regulatory region for trpE genes in E. coli.

We have adapted the simulator described in [LRSP06] and available at [RNA] to obtain sequences of terms. As an example in Figure 6 we give one (slightly shortened and simplified) terminating trajectory of the regulation process for the trpE genes (responsible for the synthesis of tryptophan) in E. coli. The regulatory region itself is presented in Figure 5.

References to the literature on RNA regulation mechanisms can be found in [LRSP06,LPRS07].

Term rewriting systems have been used in the so called Regular Model Checking framework [KMM + 01,BT02,AJMd02,ALdR05]. They have been successfully applied to the analysis of parameterized systems [BT02,AJMd02,ALdR05] and multithreaded programs [BT02,BT03,Tou05]. However, in the regular model checking framework, the rewriting rules are not probabilistic. This work constitutes the first step towards the extension of the regular model checking framework with probabilistic rewriting rules. This would allow for example the analysis of probabilistic parameterized systems and probabilistic multithreaded programs.

Rewriting systems have also been used in articles [BIK06,BCC + 03] to model chemical reactions. Compared to our work, the rewriting systems considered in [BIK06,BCC + 03] are not probabilistic. Moreover, these works consider the modeling of chemical reactions whereas we consider modeling of RNA secondary structure.

Finally, probabilistic term rewriting systems have also been considered in [BH03,BK02,KSMA03]. But in these works, the symbols are of fixed arities and the terms are ordered, whereas in our framework, the symbols have arbitrary arities and the terms are not ordered. Moreover, as far as we know, this is the first time that probabilistic term rewriting systems are used to model attenuation regulation.

Conclusions and perspectives

We have established that the framework of probabilistic term rewriting systems provides compact and structured description of detailed models of RNA regulation.

We intend to continue exploration of this framework. The most important task consists in the development of adequate data structures and algorithms, as well as approximation and abstraction methods for analysis of this kind of models. The next step would be a massive computational experimentation, the biological interpretation of results and validation of results by real biological data.

→ •(g) → •(f (e)) → •(d) → 40, 62 (c) * → •(b)→•(b)→ 40, 63 (b, h)→•(c) * → •(h)→•(c, h)→•(c)→ 40, 64 (c, h)→•(h)→•(b, h)→•(b) * → •(g)→ •(f (e))→•(b, h)→ 40, 65 (h)→•(b)→•(b, h)→ 40, 66 (h)→•(b)→•(b, i)→•(i)→•(c, i)→•(c, h)→•(c) * → 40, 67 (c, h)→ •(c, i)→•(c, j)→•(c, k)→•(b) * → •(k)→•(f (e))→•(b, h)→•(b, i)→•(b, j)→•(b, k)→•(b, i)→ 40, 68 (i)→•(b)→•(b, l)→•(b, h)→ •(b, j)→•(b, k)→•(h)→•(b, l(h))→•(k)→•(l(h))→•(l)→•(l(h))→ 40, 69 (h)→•(l)→•(b, l(h))→•(b, h)→•(b, l)→•(b)→ •(l(h))→ 40, 70 (h)→•(l)→•(b, l(h))→•(b, h)→•(b, l)→•(c, h)→•(m(h))→•(b)→•(b, i)→•(b, j)→•(b, k)→•(m) * → •(k)→ •(b, l(h))→ 40, 71 (l(h))→•(b, h)→•(b, l)→•(h)→•(l)→•(b)→•(b, l(h))→ 40, 72 (l(h))→•(b, h)→•(b, l)→•(h)→•(l)→ •(b)→•(b, i)→•(b, j)→•(b, n)→•(b, k)→•(b, o)→•(b)→ 40, 73 (b, l) * → •(b, o)→•(l)→•(b, l(h))→•(h)→•(l(h))→•(b, l(h))→ 40, 74 (l(h))→•(b, h)→•(b, l)→•(h)→•(l)→•(b)→•(b, h)→ 40, 75 (h)→•(b)→•(b, l(h))→•(l(h))→•(b, l)→•(l)→•(b, i)→ •(b, j)→•(b, n)→•(b, k)→•(b, o)→•(k)→•(c, h)→•(m(h))→•(c)→•(l(h))→ 40, 76 (h)→•(l)→•(b, l(h))→•(b, h)→•(b, l)→ •(b)→•(b, h, p)→•(h, p)→•(b, p)→•(p)→•(c, h, p)→•(c, p)→•(d, p)→•(f, p)→•(e, p)→•(g, p)→•(b, i)→•(b, j)→•(b, n)→ •(b, q)→•(b, k)→•(b, o)→•(c, h)→•(m(h))→•(b, h, p)→ 40, 77 (h, p)→•(b, p)→•(b, h)→•(p)→•(b)→•(h)→•(c, h, p)→ •(c, p)→•(d, p)→•(f, p)→•(e, p)→•(g, p)→•(b, l(h))→•(l(h))→•(b, l)→•(l)→•(b, i) * → •(b, o)→•(q)→•(l(h))→ 40, 78 (h)→ •(l)→•(b, l(h))→•(b, h)→•(b, l)→•(c, h)→•(m(h))→•(h, p)→•(m)→•(m(i))→•(m(j))→•(m(k))→•(m(h))→ 40, 79 (h)→ •(m)→•(c, h)→•(b, h)→•(l(h))→•(h, p)→•(l)→•(b, l(h))→•(b, l)→•(b)→•(b, h, p)→•(b, p)→•(p)→•(c, h, p)→•(b, h, p)→ 40, 80 (h, p)→•(b, p)→•(b, h)→•(p)→•(h)→•(c, h, p)→•(r(h, p))→•(b)→•(c, p)→•(c, h)→•(c)→•(d, p)→•(r(p))→•(f, p)→ •(e, p)→•(g, p)→•(b, l(h))→•(b, h, p)→ 40, 81 (h, p)→•(b, p)→•(b, h)→•(p)→•(b)→•(b, l) * → •(b, s)→•(h)→•(c, h, p)→ •(r(h, p)) → •(b, l(h)) → •(l(h)) → •(b, h, p) → 40, 82 (h, p) → •(b, p) → •(b, h) → •(p) → •(b) → •(c, p) → •(d, p) → •(r(p)) → •(f, p) → •(e, p) → •(g, p) → •(r) → •(r(e, p)) → •(r(g, p)) → •(r(h, p)) → •(r(e)) → •(r(h)) → •(h) → •(r(m(h))) → •(m(h)) → •(r(m))→•(m)→•(r(m(i)))→•(r(m(j)))→•(r(m(k)))→•(r(m(h)))→ 40, 83 (m(h))→•(r(h))→•(r(m))→•(r(m(h)))→ 40, 84 (m(h))→•(r(h))→•(r(m))→•(r(m(h)))→

  13, 13 () → 16, 27 (a) → 19, 27 () → 40, 50 (b) → 40, 51 (c) → •(b) * → 40, 61 (c) → •(d) *

  Figure6. A simulation result: one typical terminating trajectory for classical attenuation regulation of trpE genes in E. coli. Notations: → means one rewriting;

only on its "active" part called window, as we will see below

Acknowledgments

The authors are thankful to Sergey Pirogov, Konstantin Gorbunov and Lev Rubanov for a valuable discussion. Lev Rubanov has also provided assistance in use of the Rnamodel tool. Oleg Zverkov has helped us in preparing computer graphics for this article.