extfieldsofstudy:
Physics
Computer Science
Biology

text:
Complex Reaction Kinetics in Chemistry: A unified picture suggested by Mechanics in Physics

Complex biochemical pathways or regulatory enzyme kinetics can be reduced to chains of elementary reactions, which can be described in terms of chemical kinetics. This discipline provides a set of tools for quantifying and understanding the dialogue between reactants, whose framing into a solid and consistent mathematical description is of pivotal importance in the growing field of biotechnology. Among the elementary reactions so far extensively investigated, we recall the socalled Michaelis-Menten scheme and the Hill positive-cooperative kinetics, which apply to molecular binding and are characterized by the absence and the presence, respectively, of cooperative interactions between binding sites, giving rise to qualitative different phenomenologies. However, there is evidence of reactions displaying a more complex, and by far less understood, pattern: these follow the positive-cooperative scenario at small substrate concentration, yet negative-cooperative effects emerge and get stronger as the substrate concentration is increased. In this paper we analyze the structural analogy between the mathematical backbone of (classical) reaction kinetics in Chemistry and that of (classical) mechanics in Physics: techniques and results from the latter shall be used to infer properties on the former.

A. The Chemical Kinetics background
The mathematical models that describe reaction kinetics provide chemists and chemical engineers with tools to better understand, depict and possibly control a broad range of chemical processes (see e.g., [27,48]). These include applications to pharmacology, environmental pollution monitoring, food industry, etc. In particular, biological systems are often characterized by complex chemical pathways whose modeling is rather challenging and can not be recast in standard schemes [2,3,9,10,20,23,28,36,38,42,43,55,56] (see also [57,59,60,64] for a different perspective). In general, one tries to split such sophisticated systems into a set of elementary constituents, in mutual interaction, and for which a clear formalization is available [4,6,39,44,58,62].
In this context, one of the best consolidated, elementary scheme is given by the Michaelis-Menten law. This was originally introduced by Leonor Michaelis and Maud Menten to describe enzyme kinetics and can be applied to systems made of two reactants, say A (the binding molecule or, more generally, the binding sites of a molecule) and B (the free ligand, i.e., the substrate), which can bind (and unbind) to form the product AB. If we call S the concentration of free ligand, Y the saturation function (or fractional occupancy), namely the fraction of bound molecules (Y ∈ [0, 1]), and, accordingly, 1 − Y the fraction of the unbound molecules, under proper assumptions, one can write where k is the proportionality constant between response and occupancy (otherwise stated, it is the ratio between the dissociation and the association constants). In particular, as standard, it is assumed that the (a) the reaction is in a steady state, with the product being formed and consumed at the same rate, (b) the free ligand concentration is in large excess over that of the binding molecules in such a way that it can be considered as constant along the reaction, (c) all the binding molecules are equivalent and independent. Also, the derivation of the Michaelis-Menten law is based on the law of mass action.
Reshuffling the previous equation we get Y = S/(S + k) which allows stating that k is the concentration of free ligand at which 50% of the binding sites are occupied (that is, when S = k, then Y = 1/2). Thus, denoting with S 0 the half-saturation ligand concentration, we get This equation represents a rectangular hyperbola with horizontal asymptote corresponding to full saturation, that is Y = 1; this is the typical outcome expected for systems where no interaction between binding sites is at work [63]. This model immediately settled down as the paradigm for Chemical Kinetics, somehow similarly to the perfect gas model (where atoms, or molecules -collisions apart -do not interact) of the Kinetic Theory in the early Statistical Physics [49]. Nevertheless, deviations from this behaviour were not late to arrive: the most common phenomenon was the occurrence of a positive cooperation among the binding sites of a multi-site molecule. Actually, many polymers and proteins exhibit cooperativity, meaning that the ligand binds in a non-independent way: if, upon a ligand binding, the probability of further binding (by other ligands) is enhanced, the system is said to show positive cooperativity.
To fix ideas, let us make a practical example and let us consider the case of a well-known protein, i.e., the hemoglobin. This is responsible for oxygen transport throughout the body and it ultimately allows cellular respiration. Such features stem from hemoglobin's ability to bind (and to dislodge as needed) up to four molecules of oxygen in a nonindependent way: if one of the four sites has captured an oxygen molecule, then the probability that the remaining three sites will capture further oxygen increases, and vice versa. As a result, if the protein is in an environment rich of oxygen (e.g., in the lungs), it readily binds up to four molecules of oxygen, and, as much readily, it gets rid of them when crossing an oxygen-deficient environment. To study quantitatively its behaviour one typically measures its characteristic input-output relation. This can be achieved by considering a set of M elementary experiments where these proteins, in the same amount for each experiment, are prepared in a baker and allowed to bind oxygen, which is supplied at different concentrations S i for different experiments (e.g., S 1 < S 2 < ... < S M ). We can then construct a Cartesian plane, where on the abscissas we set the concentration of the ligand S (oxygen in this case, i.e. the input) while on the y-axes we put the fraction of protein bound sites Y (the saturation function, i.e., the output). In this way, for each experiment, once reached the chemical equilibrium, we get a saturation level and we can draw a point in the considered Cartesian plane; interpolating between all the points a sigmoidal curve will emerge (see Fig. 1). Archibald V. Hill formulated a description for the behavior of Y with respect to S: the so-called Hill equation empirically describes the fraction of molecules binding sites, occupied by the ligand, as a function of the ligand concentration [7,8,21,65]. This equation generalizes the Michaelis-Menten law (2) and reads as where n H is referred to as Hill coefficient and can be interpreted as the effective number of binding sites that are interacting each other. This number can be measured as the slope of the curve log[Y /(1−Y )] versus log(S), calculated at the half-saturation point. Of course, if n H = 1 there is no cooperation at all and each binding site acts independently of the others (and, consistently, Michaelis-Menten kinetics is restored), viceversa, if n H > 1, the reaction is said to be cooperative (just like in hemoglobin), and if n H 1 the cooperation among binding sites is so strong that the sigmoid becomes close to a step function and the kinetics is named ultra-sensitive.
The Michaelis-Menten law together with the extension by Hill, provided a good description for a bulk of chemical reactions, however, things were not perfect yet. For instance, some yeast's proteins (e.g., the Glyceraldehyde 3-Phosphate Dehydrogenase [24]) produced novel (mild) deviations from the Hill curve: for these enzymes, the cooperativity of their binding sites decreases while increasing the ligand concentration. The following work by Daniel E. Koshland allowed understanding this kind of phenomenology by further enlarging the theoretical framework through the introduction of the concept of negative cooperativity. In fact, in the previous example, beyond the positive cooperation between the binding sites there are also negative-cooperative effects underlying. Their effective action is to diminish the overall binding capabilities of the enzyme and thus to reduce the magnitude of its Hill coefficient. For the latter we report two fits: Dashed line is the result obtained by constraining the system to be cooperative but not-ultrasensitive (that is, J ≤ 1), while solid line is the best fit (without constraints) which yields to J ∼ 1.1, hence a "first order phase transition" in the language of statistical mechanics. The relative goodness of the fits are R 2 coop ∼ 0.85 and R 2 ultra ∼ 0.94, confirming an ultra-sensitive behavior. The tables in the bottom present the value of J derived from the best fit and the resulting nH ; the estimate of the Hill coefficient taken from the literature is also shown for comparison. This figure was presented in [7].

B. The Mechanics background
The progressive enlargement of a theoretical scaffold to fit the always increasing amount of evidences is a common feature in the historical development of scientific disciplines [34,41]. This is the case also for Mechanics and, as we will see, the analogy with Chemical Kinetics goes far beyond this feature. Beyond Kinematics, which describes the motion of systems without considering their mass or the forces that caused the motion, in the seventeen-th century Newton gave a sharp description of Mechanics, in the form of laws describing how masses dynamically respond when stimulated by an external force (or moment). Here, the input is the force while the output is the motion of the body. The Newtonian dynamics has been ruling for centuries and, in fact, it was so well-consolidated that scientists, among which Giuseppe L. Lagrange, William R. Hamilton and Carl G.J. Jacobi, later reformulated the entire theory in a powerful and elegant variational flavor. The theory was overall brilliant to explain the perceivable reality, but with exceptions emerging in the limit of too little or too fast. We will focus on the latter. In the Newtonian world, if an applied force is kept constant over a mass, this will constantly accelerate, eventually reaching diverging velocities. This was perfectly consistent with the general credo that the speed of light were infinite. However, this postulate broke down in 1887 when the famous experiment by Albert A. Michelson and Edward Morley proved that such a velocity is actually finite. The next years were dense of novel approaches and ideas by many scientists, as Hendrik Lorentz and Hermann Minkowski, and culminated with the special relativity by Albert Einstein in 1905. According to this theory, no mass can exist whose velocity may diverge, the limiting speed being the speed of light. The classical Hamilton-Jacobi equations and Galilean transformations left the place to the Klein-Gordon formulation and Lorentz covariances and contravariances (the natural metric being Minkowskian) [46]. Clearly, classical mechanics was still a good reference framework for the vast majority of the data collected (much like the positive cooperativity accounted for the bulk of the empirical data in the chemical counterpart), however there were rare phenomena (e.g., a muon decay in atmosphere [18]) that required a broader scaffold which, in the opportune limits, could recover the classical one.
Despite this historical connection between Chemical Kinetics and Classical Mechanics may look weird at a first glance, as we will prove, there is a formal analogy between their mathematical representations. In the next section we will summarize the main results concerning the analogy at the classical level. More sharply, the saturation plot of classical (positive cooperative) chemical kinetics (namely the input-output relation between the saturation function and the concentration of the substrate) can be derived by a minimum action principle that is the same that holds in classical mechanics, when describing a mass motion in the Hamilton-Jacobi framework. In this parallelism, the saturation function in Chemistry plays as the velocity in Physics: thus, exactly as happens in special relativity, the velocity of the mass is bounded (by definition, the saturation function can not exceed one). Indeed, we can follow this mathematical equivalence and verify that there is actually a natural broader formulation for chemical kinetics that is exactly through the Klein-Gordon setting (rather than its classical Hamilton-Jacobi counterpart) and the theory as a whole is Lorentz-invariant. Remarkably, when read with chemical glasses, this extended relativistic setting allows for the anti-cooperative corrections that Koshland revealed in the study of the yeast enzymes, resulting in a complex mixture of positive and negative cooperation among binding sites.

II. THE STANDARD MATHEMATICAL SCAFFOLD FOR CLASSICAL COOPERATIVITY
As anticipated in Sec. I A, cooperativity is a widespread phenomenon in Chemistry and its underlying mechanisms can be multiple: for example, if the adjacent binding sites of a protein can accommodate charged ions, the attraction/repulsion between the ions themselves may result in a positive/negative kinetics; in most common cases, the bonds with the substrate modify the protein conformational structure, by influencing possible further links in an allosteric way [6,52]. Whatever the origin, cooperativity in Chemistry is a typical emergent property that directly relates the microscopic description of a system at the single binding-site level, with the macroscopic properties shown by its constituent molecules, cells, and organisms, thus the use of Statistical Physics for its investigation appears quite natural [7,63]. Usually, in Statistical Physics one is provided with an (inverse) temperature β, and with an Hamiltonian (i.e., a cost function) H(σ, J, h) describing the model at the microscopic level, namely in terms of elementary variables σ i , σ j , couplings among elementary variables J ij , and external fields h i acting over these. The goal is to obtain the free energy A(β, J, h) of the model, from which the average value of the macroscopic observables can be directly derived [63].

A. Formulation of the problem: the thermodynamical free energy
In the following we summarize the minimal assumptions needed when modelling chemical kinetics from the Statistical Physics perspective; for a more extensive treatment of this kind of modelling we refer to [5][6][7][8]63], while for a rigorous explanation of the underlying equivalence between Statistical Mechanics and Analytical Mechanics we refer to the seminal works by Francesco Guerra [32], dealing with the Sherrington-Kirkpatrick model (and then deepened in e.g., [14,16,17,29]), and by Jordan G. Brankov and Valentin Zagrebnov in [31], dealing with the Husimi-Temperley model (and then deepened in e.g., [11,12,15,30]).
• Each binding site may or may not be occupied by a ligand: this allows us to code its state (empty versus full) by a Boolean variable. For the generic i th site, we will use an Ising spin σ i = ±1, where σ i = −1 represents an empty i th site, viceversa σ i = +1 means that the i th site is occupied. Clearly, if there are overall N binding sites, i ∈ (1, ..., N ).
• It is rather inconvenient, an ultimately unnecessary, to deal with the whole set σ i , ..., σ N if we are interested in the properties of large numbers of these variables (i.e., in the so called thermodynamic limit corresponding to N → ∞). If we want to distinguish between a fully empty state σ i = −1∀i ∈ (1, ..., N ) (ordered case), a fully occupied state σ i = +1∀i ∈ (1, ..., N ) (ordered case) or a completely random case where σ i = ±1 with equal probability (disordered case), it is convenient to introduce the order parameter for these variables as the magnetization M (this term stems from the original application of the Statistical Mechanics model in the context of magnetism) that reads as the arithmetic average of the spin state, namely There is a univocal relation between the magnetization in Physics and the saturation function in Chemistry, where, we recall, we denote with Y ∈ [0, 1] the fractional occupation of the binding sites. In fact, one has [7,8] Eq.(5) constitutes the first bridge between the Chemistry we aim to describe (via the saturation function Y ) and the Physics that we want to use (via the magnetization M ).
• All the binding sites interact with the ligand by the same strength. This is a standard assumption in Chemical Kinetics [8,50,65] and it means that the diffusion of the ligands is fast enough to ensure a homogeneous solution. The concentration of free ligands is mapped into a one-body contribution H 1 in the cost-function. This term encodes for the action of an external magnetic field in such a way that, if the field acting on the i th is positive, the spin will tend to align upwards (namely this direction is energetically favored), and vice versa. This homogenous mixing assumption translates into a homogeneous external field h, and the related contribution reads as Notice that h plays as a chemical potential and, consistently, it can be related to the substrate concentration as S 0 being the value of the ligand concentration at half saturation.
Eq. (7) constitutes the second bridge between the Chemistry we aim to describe (via the ligand concentration S) and the Physics that we want to use (via the magnetic field h).
• The binding sites can cooperate in a positive manner: this can be modelled by introducing a coupling between the σ variables. The simplest mathematical form is given by a two-body contribution H 2 in the cost function. This term encodes for the reciprocal interactions among binding sites and it reads as where J ≥ 0 is the interaction strength and the sum runs over all possible pairs; the normalization factor 1/N ensures the linear extensivity of the cost-function with respect to the system size. A positive value for J implies an imitative interaction among binding sites: configurations where spins tend to be aligned each others (namely where sites tend to be either all occupied or all unoccupied) are energetically more favoured and will therefore be more likely.
• Combining together the previous contributions we get the total Hamiltonian It is possible to introduce the free energy associated to such an Hamiltonian as where β is the inverse temperature in proper units. The free energy is a key observable because it corresponds to the difference between the internal energy U and the entropy S (at given temperature), i.e. A(β, J, h) = S(β, J, h) − βU (β, J, h). If we could obtain an explicit expression for A(β, J, h) in terms of the order parameter M , we could obtain an expression for the magnetization expected at equilibrium by imposing δ M A(β, J, h) = 0, in fact, this implies that we are simultaneously asking for the minimum energy and the maximum entropy.
Notice that, having stated the two bridges given by Eq.s (5) and (7), other mappings between the two fields (e.g., the relation between the coupling strength J and the Hill coefficient n H , see eq.(25) later on) emerge spontaneously as properties of the thermodynamic solutions of the problems.
B. Resolution of the problem: the mechanical action We want to find an explicit expression (in terms of M ) for the free energy defined in eq. (10). To this task let us rename −βJ = t and βh = x and let us think of these fictitious variables as a time and a space, respectively. Thus, we can write the free energy as where we also wrote i<j σ i σ j ∼ (1/2) i,j σ i σ j , which implies vanishing corrections in the thermodynamic limit. If we work out the spatial and temporal derivatives of the free energy (12) we obtain where the average · t,x for a generic observable O depending on the spin configuration is defined as and, posing t = −βJ and x = βh, the Boltzmann average for the original system (9) is recovered and this shall be simply denoted as · If we now introduce a potential V (t, x), defined as half the variance of the magnetization, i.e., we see that, by construction, the free energy of this model obeys the following Hamilton-Jacobi equation and therefore A(t, x) is also an action of Classical Mechanics. We can simplify the previous equation by noticing that, for large enough volumes, the magnetization is a self-averaging quantity [14,63], thus in the infinite volume limit the potential must vanish, that is, lim N →∞ V (t, x) = 0. Here, we are restricting to large volumes and we are therefore left with a Hamilton-Jacobi equation describing a free propagation; since the potential is zero, the Lagrangian L coupled to the motion is just the kinetic term that is the analogous of the classical formula L = 1 2 mv 2 where the mass m is set unitary (i.e., m = 1), and the role of the velocity v is played by the average magnetization M . Solving the Hamilton-Jacobi equation is then straightforward: the solution is formally written as The evaluation of the Cauchy condition A(t = 0, x = x 0 ) is trivial because, at t = 0, the coupling between variables disappears (see eq.(10)), while the integral of the Lagrangian over time reduces to the Lagrangian times time (as the potential is zero). Pasting these two contributions together we obtain Finally, noticing that the equation of motion is a Galilean trajectory as recasting the solution back in the original variables, i.e. t = −βJ and x = βh, we get the free energy associated to this general positive cooperative reaction: By extremizing A(β, J, h) with respect to M we get This result recovers the well known self-consistency equation for the order parameter of the Curie-Weiss model in Statistical Mechanics [14,63].

C. Chemical properties of the physical solution
The self-consistent equation in Eq. (22) is an input-output relation for a general system of binary elements, possibly positively interacting, under the influence of an external field: the input in the system is the external field h and the output is the magnetization M . We can rewrite Eq.(22) in a chemical jargon by using the bridges coded in the Eq.s (5) and (7) and fixing, for the sake of simplicity, S 0 = 1, that is Before proceeding, we check that if cooperation disappears (i.e., binding sites are reciprocally independent), the Michaelis-Menten scenario is recovered. Posing J = 0 in the equation above we get that is (apart a constant factor that can be re-introduced by taking S 0 = k, rather than 1) the Michaelis-Menten equation (see Eq. 2). One step forward, we now take into account the coupling J and relate it to the Hill coefficient n H . The latter is defined in Chemistry as the slope of Y (S) at half saturation (i.e., when Y = 1/2), and we can obtain its expression following this prescription by using Eq.(23), namely We note that as J → 0 we get, as expected, n H → 1: if there is no cooperation between binding sites, the Hill coefficient must be unitary; further, the stronger the coupling J and the (hyperbolically) larger the value of the Hill coefficient. In particular, for J → 1 the kinetics gets ultra-sensitive and discontinuities emerge. We remark that, with simpler statistical mechanics model as linear chains of spins, phase transitions are not allowed, hence ultra-sensitive behavior can not be taken into account: the present framework is the simplest non-trivial scheme where all these phenomena can be recovered at once (see Fig. 1 and [7] for more details on ultra-sensitive kinetics). Also, it is worth highlighting the full consistency between our treatment of ultra-sensitive kinetics and more standard ones as for instance reported in [48] (see eq.5.17 therein), where the expression for the Hill coefficient can be translated into our formulation as One we see that for M → ±1 the Hill coefficient diverges, which is the signature of an ultra-sensitive behavior: this is perfectly coherent with our approach where, in that limit, the input-output relation (see the hyperbolic tangent (23)) becomes a step function. However, as mentioned in the Introduction, this theory has its flaws, in Chemistry as well as in Mechanics. Regarding the former, the complex picture of yeast's enzymes evidenced by Koshland [24,25], where positive and negative cooperativity appear simultaneously (and with the anti-cooperativity effect getting more and more pronounced as the substrate concentration is raised), still escapes from this mathematical architecture. Further, from the mechanical point of view, two weird things happen: the velocity M is bounded by c = 1, while in Classical Mechanics the velocity may diverge; further, if we look at the Boltzmannfaktor in the free energy (Eq. 12), this reads as exp N −tM 2 /2 + xM and, recalling that the kinetic energy in this mechanical analogy reads as M 2 /2 (the mass is unitary, thus velocity and momentum coincide), we are allowed to interpret A(β, J, h) as a real action. From this perspective, the exponent can be thought of as the coupling between the stress-energy tensor and the metric tensor: a glance at the form of the Boltzmannfaktor reveals that the natural underlying metric is (−1, +1) rather than (+1, +1) as in classical Euclidean frames, or in other words, it is of the Minkowskian type. All these details point toward the generalization of the equivalence including special relativity. Plan of the next section is to follow the mechanical path and extend the classical kinetic energy including relativistic corrections and then to investigate its implications. We will see that in the broader, relativistic framework for chemical kinetics the deviations that Koshland explained adding an anti-cooperative interactions -beyond the cooperative onesat high ligand's doses are the chemical analogies of the deviation from classical mechanics at high velocities observed in special relativity.

A. Relativistic setting
The relativistic extension of the the Hamiltonian (9) is defined by an Hamiltonian of the form where M = 1 N N i σ i as usual. Next, we have to insert (26) into the free energy (10): where we already replaced t = −βJ and x = βh in order to work out their streaming, that read as where the Boltzmann averages . . . t,x are defined as (using the magnetization as a trial function) As before, the averages . . . t,x will be denoted by . . . whenever evaluated in the sense of thermodynamics (i.e. for t = −βJ and x = βh). By a direct calculation, it is straightforward to see that the expression (27) obeys the relativistic Hamilton-Jacobi equation where the symbol represents the D'Alambert operator and V (t, x) is the potential whose expression is chosen in order to make the equation valid by construction and, this time, it is automatically Lorentz invariant. If the functional A(t, x) is sufficiently smooth (that is, its derivatives are regular functions of t and x), in the thermodynamic limit, we have lim N →∞ V (t, x) = 0, hence in this high-volume limit we are left with which is the Klein-Gordon equation for a free relativistic particle with unitary mass in natural units (m 0 = 1). In relativistic mechanics, the stress energy tensor of this particle is defined as where v is the classical velocity of the particle, γ = 1/ √ 1 − v 2 and E = γm 0 = γ is the relativistic energy. In addition, the contravariant momentum is expressed through the action by the following equation Comparing (32) and (33), it is immediate to identify the magnetization as the relativistic dynamical variable: while the Lorentz factor is In the thermodynamic limit, the particle is free and its trajectories are the straight lines x = x 0 + vt. Since the relativistic Lagrangian L = −γ −1 is constant along these classical trajectories, the free energy can be computed as: Setting t = −Jβ and x = βh, we finally get an explicit expression for the free energy: Requiring that the free energy is extremal with respect to the magnetization (that from a thermodynamical perspective can be seen as the simultaneous effect of the minimum energy and the maximum entropy principles, and, from a mechanical perspective as the minimum action principle), the associated self-consistency equation becomes

B. The classical limit from a chemical perspective
Reading the self-consistency (38) in chemical terms, that is using the bridges Eq.s (5,7), we obtain We can now check whether, under suitable conditions, this broader theory recovers the classical limit. First, we notice that under the assumption of no interactions among binding sites (i.e., J = 0) and replacing h = 1 2 log(S/S 0 ), the Michaelis-Menten behaviour is recovered. This can be shown by rewriting Eq. 39 as where we also shifted S/S 0 → S for simplicity. For J = 0 the previous equation reduces to Y (S) = S/(1+S). Further, taking the classical limit, at the lowest order, we have the following expansions such that (38) reduces to (22), in the physical context, and to (23), in chemical context. Clearly, also the slope at Y = 1/2 is preserved hence, in the classical limit, we recover the expected expression for Hill coefficient (see Eq. 25), namely

C. Beyond the classical limit
To understand why we expect variations with respect to the Hill paradigm at relatively large values of the substrate concentration, we must check carefully the relativistic self-consistency (38). Let us assume we are working at nottoo high velocities (that is M < 1) and we can expand the argument inside the hyperbolic tangent, in particular, approximating 1/( The relativistic effects in chemical kinetics become transparent in this way: if we look at the field felt by the binding sites (i.e., the argument inside the hyperbolic tangent), we see that, beyond the standard Hill term βJ M (that positively pairs binding sites together) another term appears that, this time, negatively pairs binding sites together. Retaining this level of approximation, we could write an effective Hamiltonian to generate Eq. (43) that reads as hence, beyond the two-body positive coupling coded by the first term, another four-body negative coupling appears. The latter is responsible for the deviation from the classical paradigm and these deviations are in full agreement with the Koshland generalization toward the concept of mixed positive and negative cooperativity [24]. In particular, we can see at work the entire reasoning of Koshland who pointed out how, at large enough substrate concentration, the positivity of the reaction diminishes. In fact, for M ∼ 0 no relativistic effect can be noted. By increasing S (the input in the system), we get a growth in M (the output in the system): the latter raises in response of S and it is enhanced because of the two-body term in the effective Hamiltonian (44), the four-body term still being negligible. As S keeps on growing, M increases as well, up to a point where it reaches high enough values such that, from now on, also the four-body term inside the effective Hamiltonian (44) becomes relevant. At this point, a novel, anti-cooperative effect is naturally induced in the reaction and it yields to a reduction of the Hill coefficient. In the next analysis these qualitative remarks shall be addressed in more details. We focus on the definition of the Hill coefficient based on the Hill equation This equation accounts for the possibility that not all receptor sites are independent: here n H is the average number of interacting sites and the slope of the Hill plot. The latter is based on a linear transformation made by rearranging Eq. (45) and taking the logarithm: Thus, one plots log Y /(1 − Y ) versus log S, fits with a linear function and the resulting slope, calculated at the halfsaturation point, provides the Hill coefficient. As already underlined, the Michaelis-Menten theory corresponds to n H = 1 and any deviations from a slope of 1 tell us about deviation from that model. For the (classical and relativistic) models analyzed here (coded in the Hamiltonians (9) and (26)) we can estimate the slope n H directly from the self-consistency equations (23) and (39). Let us start with the classical model. We preliminary notice that Therefore, we just need to evaluate dY /d log(S) in Y = 1/2, which reads as Posing x = dY /d log(S)| Y =1/2 and noticing that S = 1 when Y = 1/2, we have By plugging this result in Eq. 47, we finally have One can see that when J = 0 the Hill coefficient is unitary as expected for non-cooperative systems, when J > 0 the coefficient is larger than 1, indicating that receptors are interacting, and when J < 0 the coefficient is smaller than 1, as expected for negative cooperativity.
Let us now move to the relativistic model. Again, we just need to evaluate dY /d log(S) in Y = 1/2, which, recalling (39), reads as Exploiting the fact that S = 1 when Y = 1/2, the previous expression simplifies as Thus, we can write Note that n class H /n rel H < 1, confirming that the relativistic correction weakens the emerging cooperativity.

D. Further robustness checks
As stressed above, for a fixed interaction coupling J, the relativistic model is expected to exhibit a lower cooperativity with respect to the classical model. In order to quantify this point we considered different quantifiers for cooperativity and we compared the outcomes for the relativistic and the classical models set at the same value of J. Let us start with the Koshland measure of cooperativity which is defined as the ratio [66]  where R 0.9 denotes the substrate concentration corresponding to a 90% saturation, while R 0.1 denotes the substrate concentration corresponding to a 10% saturation, that is, Y (R 0.9 ) = 0.9 and Y (R 0.1 ) = 0.1. In the non-cooperative case one has R 0.9 /R 0.1 = 81 and, accordingly, if the ratio is smaller than 81 (meaning that the saturation curve is relatively steep) one has positive cooperativity, while if the ratio is larger than 81 one has negative cooperativity. The advantage in using the index κ is that it can be easily measured, yet it ignores all information that can be derived from the shape of Y (S). In particular, this quantifier can be estimated starting from a Klotz plot (see e.g., Fig. 2, panel a) where the saturation function is shown versus the logarithm of the (free) ligand concentration; in the presence of positive cooperativity this plot yields to a characteristic sigmoidal curve. For the models analyzed here we can estimate R 0.9 /R 0.1 directly from the self-consistency equations (23)- (39). Starting from the classical model and, posing Y = 0.9 and Y = 0.1 we get, respectively Of course, when J = 0 we recover the value 81, when J > 0 we get R class < 81, and when J < 0 we get R class > 81. Repeating analogous calculations for the relativistic model we get and, with some algebra, log(S 0.9 ) = 2 atanh that is, Again, one can check that when J = 0 we recover the value 81, when J > 0 we get κ rel < 81, and when J < 0 we get κ rel > 81. Also, κ rel /κ class = e −16J/ √ 41+16J/5 > 1. This means that, even with this quantifier, when fixing the same coupling constant J, the emerging cooperativity is weaker for the relativistic model, as expected.
Next, let us consider the cooperativity quantifier derived from the Scatchard plot. We recall that this plot is built by showing the behavior of Y /S with respect to Y . In fact, according to the simplest scenario [67], at equilibrium, one can write where k is the proportionality constant between response and occupancy (i.e., it is the ratio between the dissociation and the association constants), and rearranging Eq. 62 we have The previous expression fits the equation of a line for Y /S versus Y , whose slope is −1/k. The advantages in using the Scatchard plot is that it is a very powerful tool for identifying deviations from the simple model, which, without deviations, is represented by a straight line. In particular, a concave-up curve may indicate the presence of negative cooperativity between binding sites, while a concave-down curve is indicative of positive cooperativity. Also, in the latter case, the maxima occurs at the fractional occupancy Y * which fulfills where σ provides another quantifier for cooperativity.
Starting from the classical model, we can build the function Y /S, by first getting S as a function of Y and can be obtained by inverting formula (23), namely By deriving Y /S with respect to Y we get which is identically equal to −1 when J = 0, monotonically decreasing with Y when J > 0 and monotonically increasing with Y when J < 0. The (possible) root therefore provides the extremal point, that is and, comparing with Eq. 64, we get We now repeat analogous calculations for the relativistic model. First, we get S as a function of Y , by inverting formula (39), namely By deriving Y /S with respect to Y we get d dY We can obtain an estimate of the value Y * corresponding to the maximum by recalling Y ≤ 1 and neglecting high-order terms. In this way we get and, comparing with Eq. 64, we get The three plots considered here (i.e., Klotz, Scatchard, and Hill), and the related estimates for the extent of cooperativity, are presented in Fig. 2. In particular, in panel d) we compare the cooperativity quantifiers for several values of J: as anticipated, in general, for a given value of J, the relativistic model gives rise to a weaker cooperativity.
We proceed our analysis by deepening the role of the coupling constant J in the binding curves related to the two models. In Fig. 3 we present the Klotz's plot (panel a), the Scatchard plot (panel b), and the Hill plot (panel c) for the relativistic and the classic models, comparing the outcomes for different values of J. As expected, the point corresponding to S = 1 and Y = 1/2 is a fixed point in each plot and, in general, the gap between the two models is enhanced when J is larger (i.e., when J is closer to 1). Also, when J is not too small, the Scatchard plot for the relativistic model displays a flex at small values of Y suggesting that, when the saturation is small, the cooperativity is not truly positive.
In the final part of this section we want to get deeper in the comparison between the classical and the relativistic models. To this aim, we solved numerically Eq. 39, for different values of S and of J, getting a set of data Y (S, J). We can think of this set of data as the result of a set of measurements where we collect the saturation value at a given substrate concentration. Now, assuming that in this experiment we have no hints about the underlying cooperative mechanisms, we may apply the formulas for the plain positive cooperativity and infer the value of J. More practically, we calculate numerically Y from the relativistic model for different values of S and of the coupling strength, referred to as J rel for clarity. Next, we manipulate the set of data Y (S, J rel ) by inverting the formula in Eq. 23: as the value of S is assumed to be known, we can derive the coupling strength, referred to as J class , expected within a classical framework. In this way, we can compare the original coupling constant J rel with the inferred one J class . We can translate this procedures in formulas as follows:  (39). From this data we inferred the expected classical coupling J class by inverting the self-consistent equation (23). We repeated the same operations for several values of S and J rel . In the left panel we show the inferred J class versus the fixed J rel : different colors represent different values of S and the identity function is also shown for reference (dashed, black curve). Notice that, in general J class < J rel and the inequality is enhanced as S grows. In the right panel we show a contour plot for the ratio J class /J rel versus h = log(S)/2 and J rel . Again, one can notice that, in general, J class /J rel < 1 and this inequality is enhanced for relatively large values of S.
In Fig. 4 (panel a) we plot J class versus J rel , for different values of S. Notice that the two parameters are related by a linear law, whose slope is smaller than 1 and decreases with S. This confirms that the relativistic model yields to a weak cooperativity. The negative contributions in the relativistic model get more effective when J rel and S are large, as further highlighted in Fig. 4 (panel b).

IV. CONCLUSIONS
The rewards in the overall bridge linking Chemical Kinetics and Analytical Mechanics are several, both theoretical and practical, as we briefly comment. The former lie in a deeper understanding of the mathematical scaffold for modelling real phenomena: it is far from trivial that the description of chemical/thermodynamical equilibrium is formally the same as the mechanical one. In particular, the self-consistency relation (38) that emerges from the thermodynamic principles (in fact, it stems from the requirement of simultaneous entropy maximization and energy minimization) also turns out to be, in the mechanical analogy, a direct consequence of the least action principles δA(t, x) = 0. This means that the stationary point corresponds to a light perturbation of the evolution of the system in the interval [0, t]. Explicitly, we shift infinitesimally M t,x → M t,x + δ M t,x , then 0 = δA(t, x) = ∂A(t, x) from which (38) is recovered (as usual by setting t = −Jβ and x = βh), since this holds for all variations δ M t,x . Even more exciting, still by the theoretical side, is the realization of the complexity of systems presenting mixed reaction (i.e., where both positive and negative cooperativity are simultaneously at work) and the possible applications in information processing, as we are going to discuss. First, let us clarify that in the Literature we speak of complex network or complex system with (mainly) two, rather distinct, meanings: in full generality, let us consider an Hamiltonian as H(σ, J) = i<j J ij σ i σ j and let us write the two-body coupling matrix as J ij = A ij W ij , where A is the adjacency matrix, accounting for the bare topology of the system (its entry A ij is 1 if there is a link connecting the related nodes (i, j), which are therefore allowed to interact each other, and it is zero otherwise) and W is the weight matrix, accounting for the sign and the magnitude of the links (i.e. the type of interactions among binding sites).
Dealing with A, networks where the topology is very heterogeneous (e.g., the distribution of the number of links stemming from a node has a power-law scaling) are called complex networks, as it is case for the Barabasi-Albert model [13]. Dealing with W , networks where the entries of the weight matrix are both positive and negative are termed complex systems, as the Sherrington-Kirkpatrick model [51] for the so called spin-glasses. Crucially, spin glasses spontaneously show very general information-processing skills and computational capabilities: for instance Hopfield neural networks [35] and restricted Boltzmann machines [1] -key tools in Artificial Intelligence (respectively in neural networks and machine learning)-are two types of spin-glasses and it is with this last definition of complexity that we now can read the information processing capabilities of the elementary reactions we studied. For a given macromolecule under consideration, we could paste each binding sites on a node and draw the links among nodes that are interacting: if two nodes are correlated (they show positive cooperativity), their relative interaction is positive, while if two nodes are anti-correlated (they show negative cooperativity), their relative interaction is negative. Dealing with mixed reactions we have to deal with spin glasses and we can thus assess how much information has been processed in a given reaction by evaluating the amount of information processed in its corresponding spin-glass representation, using our bridge. We have already started this investigation in [5][6][7]. Finally, from a practical perspective, in the classical limit (i.e., for simple reactions) we have an explicit expression that directly relates the Hill coefficient n H , which can be measured experimentally, to the interaction coupling J assumed in the model, namely n H = 1/(1 − J). This allows designing specific models and very simple validations (at least at this coarse-grained level) and gives a new computational perspective by which analyze already developed ones (see e.g. [26,33,37,40,47]. Then, regarding complex reactions, the puzzling scenario evidenced by Koshland, finally finds out a simple descriptive framework that, crucially, also recovers to the standard cooperative scenario in the proper limit: full coherence among various, apparently antithetic, results is obtained within a unique framework.
The author declare that there is no conflict of interest regarding the publication of this paper.


METHYLENETETRAHYDROFOLATE REDUCTASE C677T GENE POLYMORPHISM AND PROSTATE CANCER RISK

Objective: The single nucleotide polymorphism C677T of the methylenetetrahydrofolate reductase (MTHFR) gene encodes a thermolabile enzyme. This polymorphism was found to be implicated in cancer susceptibility. In this study, we analyzed the distribution of the MTHFR C677T polymorphism in two cohorts; patients and controls native of East of Algeria to explore the possible association between this polymorphism and prostate cancer susceptibility. Methods: Our examination has been conducted in 98 cases and 98 healthy controls. Genotyping was realized by polymerase chain reaction-restriction fragment length polymorphism method. Results: Compared with CC homozygous, the CT heterozygous was found to have a significantly increased risk of prostate cancer (p=0.04; odds ratio [OR]=2.01, 95% confidence interval [CI]: 1.02–3.95). However, no statistically significant difference was observed concerning the TT homozygous (p=0.74; OR=1.25, 95% CI: 0.51–3.04). Conclusion: Our results indicate that the genotype CT is a risk factor for prostate cancer in East of Algeria.

INTRODUCTION
Prostate cancer is one of the most common malignant cancer in terms of prevalence and the second leading cause of death among men [1][2][3]. The methylenetetrahydrofolate reductase (MTHFR) is a key enzyme in the metabolism of folates, which are crucial for deoxyribonucleic acid (DNA) methylation, nucleotide synthesis, and genomic integrity [4]. The MTHFR catalyzes an irreversible conversion of the substrate 5,10-methylenetetrahydrofolate (CH 2 -H 4 folate) to 5-methyltetrahydrofolate (CH 3 -H 4 folate) [5,6]. The 5,10-MTHF ensures the synthesis of deoxythymidylate monophosphate (dTMP) from deoxyuridylate monophosphate (dUMP). The increases of dTMP in the nucleotide provisions underrate the misincorporation of the dUMP in DNA sequence, so it decreases double-strand breaks. The 5-MTHF ensures methylation of the amino acid homocysteine to methionine, a reaction by which, the universal donor of the onecarbon group is synthesized, the S-adenosyl methionine (SAM), essential for methylation reactions [6,7]. DNA methylation profile change could cause oncogene activation and genomic instability after hypomethylation or influence on tumor suppressor genes following hypermethylation and silences their expression, which can provide prostate tumor development [6,8,9]. The variant MTHFR C677T leads to the substitution of alanine to valine, which generates a thermolabile enzyme, with approximately 35 % of the reduction in activity in heterozygous, and 70% in homozygous [10,11]. The MTHFR C677T polymorphism has been related to many cancers such as acute lymphoblastic leukemia of adult [12], the digestive system [13,14], differentiated thyroid carcinoma [15], breast cancer [16], and prostate cancer [17]. The aim of this study was to examine the relationship between the MTHFR C677T polymorphism and prostate cancer risk in the East Algerian population.

Subjects of study
We conducted a case-control study for this purpose we recruited 98 prostate cancer cases and 98 controls. All patients were histologically diagnosed with prostate cancer; thus, all stages of this tumor development have been included [18] blood samples were collected at the Uro-Nephrology Hospital "the Department of Urology and Renal Transplantation," Constantine, Algeria. Control subjects were volunteers, healthy men, without a family history of prostate cancer and they were free of cancer at the time of blood-sample. Patients and controls are native to different regions of East of Algeria. Prostate cancer case and healthy control were matched by age (±5 years). Our research has been approved by the local Ethics Committee. All participants gave their written informed consent for approximately 7 ml peripheral blood sample, after the interview.

Genotyping
Genomic DNA was extracted from leukocytes, using salting out procedure [19]. MTHFR genotyping was performed for the single nucleotide polymorphism (SNP) C677T by polymerase chain reaction (PCR)-restriction fragment length polymorphism (RFLP). The primer sequences used were Forward: 5'-TGAAGGAGAAGGTGTCTGCGGGA-3' and Reverse: 5'-AGGACGGTGCGGTGAGAGTG-3'. PCR conditions were: Initial denaturation cycle at 94 °C for 5 min, 30 cycles, each cycle contains three stages: Denaturation at 94°C for 30 s, hybridization at 65°C for 30 s, and elongation at 72°C for 40 s, and finally a terminal elongation at 72°C for 10 min. 30 µl of PCR product was digested with restriction enzyme HinfI (Promega, Madison, WI) at 37°C overnight, and digested products were separated on 3% Agarose gel. After that, fragments were visualized using ultra-violet light.

Statistical analysis
Using Epi Info Version 6, differences in allelic and genotypic frequencies have been tested. *p˂0.05 was considered statistically significant.

RESULTS
In our study, patients and controls were aged between 49 and 89 years, with an average age of 70.33±7.65 years. General characteristics of the study population represented by; the age of the two cohorts distributed every 10 years, tumor extension, Prostate Specific Antigen (PSA) levels at the time of diagnosis and the existence of a family history of prostate cancer have been shown in Table 1.
Almost all patients have been diagnosed at age 60 and more. Also, 62.24 % of them at late stage and distant metastasis. Furthermore, 58.16 % of patients with PSA levels above 50 ng/ml, underlining delayed diagnosis. Distribution of allelic and genotypic frequencies for the polymorphism MTHFR C677T between prostate cancer patients and controls and their relation with risk of prostatic carcinogenesis is presented in Table 2.
Our results showed that the 677 CT genotype was frequent between cases (p=0.04), (odds ratio [OR]: 2.01, confidence interval [CI]: 1.02-3.95), suggesting that they may be a risk factor for prostate cancer. On the other hand, for the 677 TT genotype, any statistically significant difference between cases and controls was observed (p=0.74), (OR: 1.25, CI: 0.51-3.04). Concerning the T allele frequency of the C677T SNP was 0.41 for patients and 0.36 for controls.

DISCUSSION
This paper reports for the first time in Algeria the relation between the MTHFR C677T polymorphism and susceptibility to developing prostatic cancer. Our results showed that there was an association between the MTHFR 677CT heterozygous and prostate cancer. In agreement with our results, Van Guelpen et al. [7], Marchal et al. [20], López-Cortés et al. [17], and Abedinzadeh et al. [21] reported a positive association between the polymorphism C677T and risk of prostatic cancer in the Swedish population, Spanish population, Ecuadorian, and Asians, respectively. The significant association between the CT genotype and risk of prostate cancer may be explained by the possibility to keep a sufficient stock of methionine to support neoplastic clone progression compared with the TT homozygous [20].
However, various investigations have not found any relation between this polymorphism and prostate cancer in Caucasians population [4], the Swedish population [6], and Americans [22]. As well, the metaanalysis managed by Collin et al. [23], Zhang et al. [24], and recently, Abedinzadeh et al. [21] has concluded that the C677T polymorphism, in general, has no effect on the occurrence of prostate cancer.
Besides, other authors observed that the T allele exerts a protective impact on prostate cancer risk. Like they indicated Marchal et al. [20], and Guo et al. [25] that the genotype 677 TT is a protective factor in Spanish and Asians, respectively. As well, Cai et al. [26] showed that the TT genotype has been associated with a reduced risk of prostate cancer in the Chinese population, and the T allele exerts probably a protective effect on the risk of this cancer type. Furthermore, Li and Xu [27], who had carried out a meta-analysis, have reported the same result for the T allele. Safarinejad et al. [28] have studied the association between tumoral aggressiveness and TT, CT, and CC genotypes: The comparison showed that TT homozygous reduced by more than 50% the risk of high-grade prostate cancer (Gleason score >7).
In addition, Küçükhüseyin et al. [29] proposed that the CT genotype and the T allele might be associated with decreased risk of prostate cancer among the Turks. Singal et al. [30] suggested that the CT genotype was associated with a reduced risk of prostate cancer.
To explain the protective effect of the T allele on prostate cancer, it has been proposed that the mutant 677 TT genotype was associated with significantly decreased DNA methylation status. The reduced activity of MTHFR enzyme, coded by the T allele, decreases the SAM synthesis, the headmaster donor of a methyl group for DNA methylation reactions, favor thereby tumor suppressor genes expression [20]. Second, the thermolabile enzyme increases the pool of the 5,10-methylenetetrahydrofolate, requisite for purine and deoxythymidylate triphosphate synthesis, which is the leader nucleotide, required for DNA reparation, consequently, a decrease of uracil incorporation inside DNA, reduce thereby chromosome breaks [20,31]. Furthermore, in vitro experiences approve that inhibition of MTHFR enzyme conduct to decreasing tumor development due to the limited methionine provision, which leads to stop the cell cycle on steps S and G2, damage thus their proliferative potential [20].
The contradictory results from studies of different populations may be explained by the method by which controls have been selected, the ethnic origin of cases and controls, geographical region, lifestyle and exposition to some factors [32,33].

CONCLUSION
The results suggest for the first time, but in a small sample of subjects, that the 677 CT heterozygous is a risk factor for prostatic carcinogenesis. Further researches should be undertaken that includes other possible genes that participate in the homocysteine metabolism, such as  methionine synthase. The influence of the environment, the gene-gene as well as gene-environment interactions should also be performed to clarify their possible roles in prostatic carcinogenesis.

ACKNOWLEDGMENTS
We acknowledge all patients and controls for their participation in this study. We also thank all personnel at the Department of Urology and Renal Transplantation, and personnel at the Laboratory of Biology and Molecular Genetics, for their helping.

AUTHORS CONTRIBUTION
Tellouche-Bouhouhou Samah performed the genotyping, drafted the manuscript, Mrs. Chellat-Rezgoune Djalila "supervisor" and Satta Dalila corrected and revised the manuscript, Abadi Noureddine director of Biology and Molecular Genetics laboratory, Dahdouh Abderrezak thesis advisor. The final manuscript has been approved by all authors.


extfieldsofstudy: Medicine Biology

s2fieldsofstudy: Biology Chemistry

text:
Pretreatment with probiotic Bifico ameliorates colitis‐associated cancer in mice: Transcriptome and gut flora profiling

Individuals with inflammatory bowel disease are at high risk of developing colitis‐associated cancer (CAC). Strategies to block the process from inflammatory bowel disease to CAC should be considered. In the experiment, we aim to explore the chemopreventive efficacy of the probiotic cocktail Bifico and its potential mechanism in azoxymethane and dextran sodium sulphate‐induced CAC in mice. Oral pretreatment of Bifico was adopted to evaluate its protective effect. The colorectums of 35 C57BL/6 mice were collected and examined for the degree of inflammation and tumorigenesis. Comparative 16S rRNA sequencing was carried out to observe Bifico‐target alterations in gene expression and microbiota structure. We found that pretreatment of Bifico alleviated intestinal inflammation and reduced tumor formation. Furthermore, we identified a subset of genes as potential targets of Bifico treatment, including CXCL1,CXCL2,CXCL3, and CXCL5, which are all ligands of C‐X‐C motif receptor 2 (CXCR2). The 16S rRNA sequencing showed that Bifico decreased the abundance of genera Desulfovibrio, Mucispirillum, and Odoribacter, and a bloom of genus Lactobacillus was detected. Notably, we found that an abundance of these Bifico‐target taxa was significantly associated with the expression of CXCR2 ligand genes. Our studies indicate that Bifico, given orally, can ameliorate CAC in mice through intervening with the possible link between Desulfovibrio, Mucispirillum, Odoribacter, Lactobacillus, and CXCR2 signaling.

| INTRODUCTION
It has been largely accepted that chronic inflammation is a pivotal causal component of somatic tumorigenesis and progression. 1,2 A representative example of this strong association is colonic carcinogenesis arising in individuals with IBD. 1,[3][4][5] Colorectal cancer has been well studied and it can be divided into 2 types: colitis-associated CRC and sporadic CRC. Prolonged IBD is one of the major risks for CAC initiation, characterized by CAC incidence ranging from 2% after 10 years of IBD history up to 18% after 30 years. 6,7 Considering the long course from IBD to CAC, there are opportunities for interventions to halt or even reverse the process of carcinogenesis in this high-risk population. 8 Evidence is accumulating that dysbiosis of intestinal microbiota may greatly influence the pathogenesis of IBD and CRC in both animal models and humans. [9][10][11][12][13][14] For example, bacterial phyla responsible for producing short-chain fatty acids are observed with decreased abundance in patients with IBD, which may serve as a cause or a consequence of permanent intestinal inflammation response in IBD. 9,15,16 A recent study also reported that Fusobacterium nucleatum contributes to chemoresistance in CRC through modulating the network of Toll-like receptors, microRNAs, and autophagy. 17 Thus, strategies of manipulating the gut flora to restore balance of the microbiota community may be therapeutically feasible in patients with IBD and CAC.
It has been reported by numerous studies that oral administration of probiotics VSL#3, a commercially available probiotic cocktail that contains 4 strains of genus Lactobacillus, 3 strains of genus Bifidobacterium, and 1 strain of genus Streptococcus, can reduce chronic inflammation and delay development of carcinoma in mouse models of CAC. 13,[18][19][20] However, conclusions on whether there is a therapeutic effect of probiotics on CAC are not exactly the same. The overwhelming majority of studies support the view that probiotics prevent inflammation and carcinogenesis. 18,19,21 In contrast, a recent study reported that concurrent treatment with VSL#3 enhanced tumorigenesis in an AOM/Il10 À/À mouse model. 22 Furthermore, the specific mechanisms underlying the beneficial effects of probiotics on CAC have not been fully investigated.
Bifico capsules were approved as an over-the-counter drug by the Chinese regulatory authority, the State Food and Drug Administration, in 2002. 23 Previous research has shown that Bifico could relieve gastritis induced by Helicobacter pylori and ameliorate experimental colitis caused by DSS in mice. 24,25 However, the effect of Bifico on CAC has not been studied yet.
The goal of the present study is to evaluate the protective efficacy of Bifico capsules against AOM/DSS-induced CAC in mice. In addition, we identified several significantly altered gut bacterial groups and genes that were involved in the protective effect of Bifico treatment. Our results confirm that the probiotic cocktail Bifico has the capacity to restrain colonic carcinogenesis in mice, providing evidence for the clinical application of Bifico to lower the morbidity of CAC in patients with recurrent IBD.

| Experimental procedure
The experimental timeline of this experiment is shown in Figure 1.
Thirty-five male C57BL/6 mice were randomly divided into 3 groups: the Ctrl group (n = 7), Model group (AOM/DSS only, n = 14), and Bifico group (n = 14). Modeling of CAC was carried out according to protocol. 26 The mice were injected i.p. with a single dose of AOM (10 mg/kg) on the first day and 3 cycles of DSS treatment was adopted thereafter. For each cycle of DSS, the mice in the Model and Bifico groups received water containing 2% DSS (w/v) for 7 days in a row, followed by 2 weeks of sterile water. The mice in the Bifico group were given Bifico capsules (4.2 g/kg, dissolved in 200 lL physiological saline, at least 1.

| Histopathological evaluation and immunohistochemistry
Colon biopsies were macroscopically examined for tumorigenesis and were made as "Swiss-rolls". 26 The tissues were then fixed in SONG ET AL.

| Analysis of microbiota sequencing data
Operational taxonomic units were identified with a standard of at least 97% similarity in the 16S rRNA sequence using USEARCH and chimeric sequences were identified using UCHIME. 27,28 Classification of representative reads for each OTU was carried out utilizing the RDP Classifier algorithm against the Silva 16S rRNA database with a confidence threshold of 0.7. 29,30 The unweighted pair-group method with arithmetic mean was used for hierarchical clustering, and principal co-ordinate analyses were performed with QIIME and R package software. A Venn diagram, bacterial taxonomic distributions, and Spearman's correlation coefficient were visualized using R software.
Analysis of similarities was carried out using QIIME. Linear discriminant analysis effect size was undertaken to receive representative bacterial groups of the Bifico and Model groups. 31

| Real-time qPCR
Total mRNA was extracted using TRIzol reagent (Invitrogen) and cDNA was synthesized using a PrimeScript RT reagent Kit (Takara Biotechnology, Dalian, China). Real-time qPCR was carried out using Hieff qPCR SYBR Green Master Mix (Yeasen, Shanghai, China) with the 7500 Real-Time PCR System (Applied Biosystems). The results were analyzed using the 2 ÀDDCT method. The primer sequences used in this experiment are shown in Table S1.

| Enzyme-linked immunosorbent assay
Levels of PGE 2 in the serum were assessed by double antibody sandwich ELISA Kits (Elabscience, Wuhan, China) according to the manufacturer's instructions.

| Statistics
Differences in expression of pro-inflammatory genes, tumor number, tumor average size, colon length, and richness of bacteria were analyzed by Student's t test using GraphPad Prism 7.0 software. Results were considered significant if two-tailed P-values were < .05.

| Bifico treatment reduces colonic inflammation
The chemopreventive effect of Bifico against colonic inflammation was evaluated. The body weight of mice decreased substantially F I G U R E 1 Experimental timeline for azoxymethane (AOM) and dextran sodium sulphate (DSS) mouse models of colitisassociated cancer with interventional administration of the probiotic cocktail Bifico. Ctrl, control; d, day following each periodic exposure to DSS and the mice partially retrieved their body weight when DSS treatment was suspended.
After the second exposure cycle to DSS, the body weight loss percentage of the mice in the Bifico group was significantly less extensive than the mice in the Model group (

| Transcriptional changes induced by AOM/DSS and probiotic mixture treatment
To fully explore the underlying mechanism of the Bifico effect, global transcriptional changes of mice in the Model and Bifico groups were identified by cDNA microarray analysis of colon tissues. Two-dimensional hierarchical clustering showed distinct gene expression patterns at different time points among different groups ( Figure S1).
We compared gene expression levels of the mice in the Model group with those in the Bifico group. We identified 300 differentially expressed genes between mice in the Model and Bifico groups at week 9 (experiencing 3 DSS cycles), finding that 166 genes were upregulated and 134 genes were downregulated (log 2 |fold change| ≥ 1.0; P < .05, marked as profile BIFICO). In addition, we were also interested in gene profile alterations during the process of CAC progression. We compared the levels of gene expression at week 9 to those at week 3, and found that 1155 genes were upregulated, whereas 453 genes were downregulated (log 2 |fold change| ≥ 1.0; P < .05, marked as profile CAC). We presented the above results with volcano plots (Figure S2A,B).
We then carried out biological pathway and gene ontology enrichment analyses using the Database for Annotation, Visualization and Integrated Discovery (https://david.ncifcrf.gov/) to identify possible biological pathways and biological processes associated with Bifico treatment. 32 We undertook the enrichment analysis on the gene expression of profile BIFICO, and identified nine significantly canonical pathways. Figure 4A shows the significant pathways. We then continued to investigate significant canonical pathways associ- Among the potential regulators, we selected 5 genes to validate their mRNA levels using RT-qPCR, including CXCL1, CXCL2, CXCL3, CXCL5, and CCL7. As presented in Figure 4D, fold changes of these genes in mice of the Bifico group were identified as significantly downregulated when compared with those in the Model group.

| Bifico alters composition of the mucosaadherent microbiota community
To determine the changes in colonic microbiota community, we assessed the mucosa-adherent microbiota from colon tissues by using high-throughput sequencing. Microbial community diversity was estimated using the Chao1 and Shannon indices, and the richness was assessed based on the number of OTUs. Diversity of the gut flora was not significantly altered by Bifico treatment, although a lower trend in diversity in the intestinal microbiota of Bifico-fed mice was observed ( Figure 5A). Similarity of the colonic flora among samples was revealed using a phylogenic tree. Each branch of the tree represented the microbiota structure of 1 sample, and the F I G U R E 4 Results of cDNA microarray and validation of the expression of selected genes in azoxymethane (AOM) and dextran sodium sulphate (DSS)-induced mouse models of colitis-associated cancer. A, Significant canonical pathways of genes expressed differentially in response to treatment with the probiotic cocktail Bifico. B, Significant canonical pathways in progression of colitis-associated cancer. C, Bifico target genes (n = 3). D, Validation of the expression levels of selected genes; n = 4 mice/treatment group. Data are expressed as mean AE SEM (n = 4). *P < .05, **P < .01, ***P < .001, compared to the Bifico group or control (Ctrl) group. Bifico 9, mice treated with AOM/DSS and Bifico, killed at 9 wk; Model 3, mice treated with AOM/DSS without Bifico, killed at 3 wk; Model 9, mice treated with AOM/DSS without Bifico, killed at 9 wk SONG ET AL. | 671 samples in the same group shared a high similarity ( Figure 5B). Next, we examined the shared and unique bacterial OTUs among different groups. The shared OTUs from the 24 samples are shown with a Venn diagram ( Figure 5C).
To obtain in-depth knowledge of changes of microbiota community structures caused by Bifico treatment, a principal co-ordinate analysis was carried out on Bray-Curtis dissimilarity and UniFrac weighted distance. Significant separation was revealed between Model and Bifico groups (analysis of similarity based on Bray-Curtis and weighted UniFrac distance; P-value, .006 and .014; Figure 5D,E, respectively).

| Bifico changes the abundance of specific taxa of the mucosal microbiota
Relative abundances of taxa (contributing >1%) are depicted using pie charts in Figure 6A,B. Distribution of colonic microbiota in the Ctrl, Model 3, and Bifico 3 groups is shown with bar charts in Figure S3. To explore the association between gut microbiota and the effect of Bifico, we further analyzed the changes on abundance of specific taxa. Bifico-treated mice showed a significant contraction of Deferribacteres in the mucosal microbiota ( Figure 6A, Table S2). At the genus level, Bifico treatment resulted in a significant shrinkage of Mucispirillum, Odoribacter, and Desulfovibrio ( Figure 6B, Table S3).
Notably, Bifico treatment contributed to a significant increment of Lactobacillus ( Figure 6B). We also attempted to identify the altered bacterial groups during the progression of CAC. Mucispirillum, Odoribacter, and Desulfovibrio were of increased abundance in Model 9 group, whereas Ruminiclostridium showed the opposite change ( Figure S4).
To further determine the core bacterial group on the intestinal mucosa of mice, LDA coupled with effect size measurements was carried out (P < .05, LDA score >2.0). We found that 31 bacterial groups were enriched in the Model group, and 15 taxa in the Bifico treatment group (Figure 7A,B). Of all these bacterial groups, Desulfovibrio, Mucispirillum, and Odoribacter were the core genera F I G U R E 5 Bifico probiotic alters the community composition of mucosa-adherent microbiota in models of colitis-associated cancer. Mice were treated with azoxymethane and dextran sodium sulphate and killed at 3 wk or 9 wk (Model 3 and Model 9 groups, respectively), or also received Bifico (Bifico 3 and Bifico 9 groups, respectively); n = 5 mice/group. A, Bifico did not alter the Shannon index or Chao1 index. B, Hierarchical clustering tree for each sample at the operational taxonomic unit (OTU) level. C, Venn diagram of shared and unique OTUs among different groups. The histogram represents the number of OTUs of each group. D, E, Principal coordinates analysis based on Bray-Curtis dissimilarity (D) and weighted UniFrac dissimilarity (E). Control (Ctrl) group, n = 4. n.s., not significant in the Model 9 group, whereas Lactobacillus was remarkably enriched in the Bifico 9 group (LDA score >4.0; Figure 7B).

| Changes in gene expression profiles are correlated with microbiota alterations
We selected 5 differentially expressed genes, concluded from results of the cDNA microarray, and quantified their expression changes ( Figure 4D). More importantly, we found that some Bifico-influenced taxa were associated with these genes. Correlation between the relative abundance of the bacterial genus and the mRNA fold-change of selected genes of mice among the Bifico 9 and Model 9 groups are depicted utilizing a heat map of Spearman's correlation coefficients ( Figure 7C). Taxa that were significantly associated with any of the 5 genes included Desulfovibrio, Odoribacter, Mucispirillum, and Lactobacillus ( Figure 7C, Table S4). The bacterial groups that were of decreased abundance in the Bifico group, Desulfovibrio and Odoribacter, were positively correlated with expression of genes CXCL1, CXCL2, CXCL3, CXCL5, and CCL7, and Mucispirillum was positively associated with CXCL1, CXCL3, CXCL5, and CCL7 (Table S4). In contrast, Lactobacillus was negatively associated with the expression of CXCL2 and CCL7 (Table S4). Other correlated bacterial groups are presented in Table S4. Similar correlations were observed at the OTU level (Table S5).

| DISCUSSION
It is universally accepted that alteration in the community structure of the microbiota acts as a driving force during the stepwise progression from inflammation to dysplasia to adenocarcinoma. Thus, the possibility of retraining tumor progression by manipulating the intestinal microbiota seems logical. In this experiment, we confirmed that prophylactic administration of Bifico could alleviate weight loss, reduce the multiplicity and average size of tumors, and lower the expression of pro-inflammatory genes in an AOM/DSS mouse model. These results are in agreement with numerous previous reports focusing on probiotics. [18][19][20]33 Notably, compared with other studies focusing on VSL#3, we used a lower dose of probiotic (VSL#3 vs Bifico, 1 9 10 9 c.f.u. vs 1.2 9 10 7 c.f.u./mouse each day) to achieve a similar inhibitory effect on inflammation and tumorigenesis in an AOM/DSS mouse model, which reflects that Lactobacillus, Bifidobacteria, and Enterococcus could be an alternative combination, and a high dosage of probiotics may not be indispensable.
F I G U R E 6 Bifico probiotic changes the structure of gut flora in models of colitis-associated cancer. Mice were treated with azoxymethane and dextran sodium sulphate and killed at 9 wk (Model 9 group), or also received Bifico (Bifico 9 group); n = 5 mice/group. Pie charts depict median standardized abundance (% of total) of the indicated taxa. A, B, Phylum level (A) and genus level (B) distribution of gut flora in Bifico 9 and Model 9 groups. Standardized abundance of taxa contributing >1% were included, and pair-wise comparisons were carried out using the Wilcoxon rank sum test. *P < .05 SONG ET AL.

| 673
To investigate the underlying mechanism of probiotics' preventive effect, we detected alterations in the gene expression profile and microbial community of the intestine. We compared gene expression profiles in the inflammatory phase (1 cycle of DSS) with that in the neoplastic phase (3 cycles of DSS), and canonical pathway analysis of these genes highlighted "cytokine-cytokine receptor interaction" as the most significantly influenced pathway induced by AOM/DSS, which is consistent with a previous study. 34 Most noteworthy was that this pathway also showed significance in a comparison of gene expression profiles between the Bifico group and the Model group, indicating that Bifico might exert preventive effects against tumorigenesis through reversing the expression of these genes.
CXCL1, CXCL2, CXCL3, and CXCL5 are all CXC chemokines with an Glu-Leu-Arg motif, and bind to CXCR2 to play their role. 35

CXCL1
and CXCL5 are reported to be significantly upregulated in human colorectal cancer specimens compared with their tumor adjacent tissues. 36 CXCL1 contributes to tumor angiogenesis, cell invasion, and migration in several types of cancer, including colon cancer. 5,37 In keeping with CXCL1 and CXCL5, CCL7 is also reported to enhance colon cancer progression and metastasis through epithelial-mesenchymal transition. 38 A previous study reported that deletion of CXCR2 can attenuate AOM/DSS-induced colitis-associated tumorigenesis through diminishing infiltration of granulocytic myeloid-derived suppressor cells. 35 Additionally, CXCR2 signaling can also enhance blood vessel formation in solid tumors through recruiting tumor-associated neutrophils. 39,40 In our study, the gene expression of CXCR2 showed no difference between the Model and Bifico groups. However, CXCR2 ligands CXCL1, CXCL2, CXCL3, and CXCL5, which were significantly induced by AOM/DSS, were all markedly downregulated by Bifico treatment, indicating that Bifico might interfere with the migration of numerous types of tumor-infiltrating inflammatory cells and therefore reshape the tumor microenvironment to inhibit tumor progression.
F I G U R E 7 Cladogram, linear discriminant analysis (LDA) distribution, and heat map of Spearman's correlation coefficients of gut flora in models of colitis-associated cancer. Mice were treated with azoxymethane and dextran sodium sulphate and killed at 9 wk (Model 9 group), or also received the probiotic Bifico (Bifico 9 group). A, Taxonomic representation of taxa of the mice in Model and Bifico groups. B, Histogram of the LDA scores for differentially abundant microbiota. The threshold of LDA score was 2.0 (n = 5). C, Heat map of Spearman's correlation coefficients between expression of C-X-C motif receptor 2 signaling genes and relative abundance of Bifico target taxa. *Bacterial groups with LDA score >4.0 (n = 4) Alternatively, Bifico could distinctly change the composition of the gut microbiota community attached to the colonic mucosa. As previous studies described, Desulfovibrio had the capacity to degrade sulfate and produce high levels of hydrogen sulfide, which is regarded as a genotoxic agent that might play a role in colitis and colonic tumorigenesis. 41 Desulfovibrio also displays a flourish of relative abundance in IBD colon tissues and in patients with adenoma and rectal cancer. 14,[42][43][44] However, most studies have shown that the abundance of Desulfovibrio in the intestinal tract of colon cancer patients does not increase. 45,46 Considering its role in triggering colonic inflammation, Desulfovibrio may be of larger significance in CAC initiation and progression than in sporadic colorectal cancer. Another putative mucin degrader, Mucispirillum, was detected at decreased levels in the Bifico-treated group when compared with the Model group. [47][48][49] Mucispirillum, with increasing abundance in DSS-treated mice and natural killer T cell-deficient mice, is recognized as opportunistic pathogen that may fuel inflammation. 14,[47][48][49] Odoribacter is also reported to markedly bloom in AOM/ DSS-induced CRC mice and patients with rectal cancer. 42,50 As shown in our results, the significant enrichment of Desulfovibrio, Mucispirillum, and Odoribacter in the Model group was partly removed by treatment with Bifico, suggesting that Bifico could relieve inflammation and decrease tumor formation by inhibiting these colitogenic bacteria. We also detected a tremendous expansion of the genus Lactobacillus, a member of the Bifico bacterial group. As reported by numerous studies, treatment with various strains of Lactobacillus reduces the severity of colitis and colitis-associated carcinogenesis in both mice and humans 21,51-54 Therefore, we assume that Lactobacillus is the crucial component in Bifico's chemoprevention of colitis-associated tumorigenesis. Nevertheless, we did not detect any significant alteration in abundance of Bifidobacteria or Enterococcus, 2 other components in the Bifico cocktail, probably because we tested mucosa tissues rather than stool specimens, and the vast majority of the bacteria was excreted with feces.
Interestingly, we found that these CAC-associated bacterial taxa were correlated with the expression of CXCR2 signaling genes, which indicates a mechanistic link between CXCR2 signaling gene expression and alterations of certain bacterial groups. Such an association sustains the view that the microbiota could regulate the development of tumors through shaping the immune system. 5,10 Presumably, Desulfovibrio and Mucispirillum may putatively degrade mucin, a component of the viscous-elastic mucus, and therefore impair its protection against external damage. Such a disruption of the intestinal barrier may directly lead to the consequence of increased exposure of the intestinal epithelial cells to the luminal microbiota. It is noteworthy that Desulfovibrio is reportedly capable of invading and surviving in KB cells, a kind of oral epidermoid carcinoma cell, thereby giving rise to increased production of interleukin-6 and interleukin-8. 55 We can speculate that Desulfovibrio might also trigger the inflammatory response through such a direct way in the process of CAC initiation and progression.
Notably, pathogenic bacteria groups are reported to induce expression of the inflammatory enzyme COX-2 in inflamed colonic mucosa. 56,57 Both COX-2 and COX-2-derived PGE 2 have been proven to upregulate CXCR2 ligands in intestinal mucosa and tumors. 35 Considering that the expression of COX-2 in the colon tissues and levels of PGE 2 in the serum were significantly downregulated by treatment with Bifico, we predict that COX-2 could be one of the mediators that links the altered colonic microbiota and the CXCR2 signaling axis in the Bifico intervention process.
Our results strongly suggest that the probiotic mixture Bifico exerted an efficacious chemopreventive impact on inflammationassociated CRC. Although we identified the noticeable changes in the global transcriptome and microbial community associated with probiotic treatment, the causal relationship between these changes and attenuated tumorigenesis cannot be determined. Further detailed studies are required to investigate this deeply. Most of all, our results imply a mechanistic link between the possible Bificoresponsive pathogenic taxa (Desulfovibrio, Mucispirillum, and Odoribacter) and the CXCR2 signaling axis during CAC progression, which could explain, in part, the pattern of the preventive effect of probiotic treatment.

ACKNOWLEDG MENTS
We thank all members of the Lab of Zhaoyuan Hou (Shanghai Key Laboratory for Tumor Microenvironment and Inflammation, Department of Biochemistry and Molecular Cellular Biology, Shanghai Jiaotong University School of Medicine, Shanghai) who helped with discussions and reviews of the manuscript. We thank Shanghai Sinepharm for its kind gift of Bifico capsules.

CONFLI CT OF INTEREST
The authors have no conflict of interest.


s2fieldsofstudy: Biology Medicine

text:
Molecular analysis of gastric cancer identifies genomic markers of drug sensitivity in Asian gastric cancer

Background: Gastric cancer (GC) is one of the leading causes of lethal malignancies worldwide, especially in Eastern Asia. Clinical responses to antitumor therapies are often limited to a subset of patients. Methods: To uncover new biomarkers of sensitivity and resistance to cancer therapeutics, we performed ultra-deep targeted sequencing in a cohort with 72 patients (41 with chemotherapy sensitivity and 31 with chemotherapy resistance). Results: We found that sixteen mutated cancer genes were associated with widely used agent in chemotherapy of gastric cancer. Genes identified in these study are mainly involved in activation and inactivation of cancer chemotherapeutic agents, changes of apoptosis and proliferation, drug efflux, DNA damage repair, and the tumor microenvironment. Discussion: A novel group of chemo-sensitivity related genes provided new therapeutic strategies to overcome the development and evolution of resistance to cancer chemotherapy.

Introduction
Gastric carcinoma is one of the most common cancers and the second leading cause of cancer-related death in the world [1]. Despite a decline in morbidity and mortality, the burden remains high [2]. In gastric cancer patients, the treatment outcomes depend on the stage of the tumor at presentation and the condition of the patients. Surgical resection is the only possible curative treatment for gastric cancer. Moreover, more than two-thirds of patients have metastatic or surgically unresectable disease [3]. Therefore, Chemotherapy is used to relieve the symptoms of unresectable tumors and reduce the risk of recurrence and metastasis after surgery [4]. Although chemotherapy plays an important role in the treatment of recurrent and metastatic gastric cancer, multiple drug resistance (MDR) to chemotherapy is known as the major cause of treatment failure for gastric carcinomas [5].
Recently, numerous studies have reported the mechanisms of drug resistance in gastric carcinomas [6,7]. However, the mechanisms of drug resistance in gastric cancer have not been fully elucidated.
To find mutations predicting chemotherapy benefit in gastric carcinomas, we performed ultra-deep targeted sequencing in a cohort with 72 patients (41 with a good response and 31 with a poor Ivyspring International Publisher response). We sequenced 1,574 genes involved in tumor progression and development. The purpose of this study was twofold, firstly to find therapeutic response predictive biomarkers, secondly to identify new targets in gastric carcinomas treatment.

Material and Methods
Patients and Tissue Samples 72 pairs of tumor and normal samples were collected from patients diagnosed with gastric carcinomas; all patients had received radical surgery (total or partial gastrectomy). After surgery, all patients were received chemotherapy (fluorouracil, platinum or paclitaxel). All patients were Chinese. The original clinical and survival data of patients were provided in Table 1. The ethical committee of Ren Ji Hospital, Shanghai Jiao Tong University School of Medicine approved the studies and all patients signed informed consent before recruitment to the study.
Tumor cell purity was assessed in hematoxylin and eosin (H&E) sections. At least 5 slices of 10 μm of thickness were cut from the paraffin block and tumor regions were scraped according to the assessment of tumor enriched area.

Ultra-deep targeted sequencing and reads alignment
Genomic DNA was extracted using the GenReadTM DNA FFPE Kit (Qiagen; Germany) from formalin-fixed paraffin-embedded tissues, as well as matched normal DNA, and was fragmented and hybridized according to the manufacturer's protocol. For ultra-deep targeted sequencing, we used SeqCap EZ Choice kit (Roche). Raw fastq files were checked quality control using FastQC, and were removed low quality reads to clean reads. Burrows Wheller aligner (BWA) was used to align the sequencing reads onto the hg19 reference genome with default options [8]. Samtools was used to convert the SAM file into BAM; Picard was used to remove duplicate reads [9].

Detection of somatic mutations and indels
Mutect was used to detect somatic mutations, single-nucleotide variants (SNVs), small insertion and deletion (INDELs) [10]. To obtain a list of confident variants, we selected variants whose coverage was more than 20X. Vcf2maf was used to convert a vcf file into a maf file and annotate the somatic mutations [11].

Statistical Analyses
Statistical analyses to test for correlation between mutations and clinicohistological variables such as age, sex and lymph node metastasis were performed using R/Bioconductor packages. We used the Fisher's exact test for difference between groups [13]. Using the clusterProfiler package, KEGG pathway enrichment analyses were performed with the significantly genes. The P value of 0.05 or less was considered statistically significant.

Tissue Samples and Clinicopathological Data
Seventy-two cases of gastric carcinomas patients (forty-one chemo-sensitive and thirty-one chemo-resistant) were analyzed (Table 1). In our cohort, the mean age was 61 years (range 33-80; SD 9.4), 57 patients were male (79%). With a median follow up of 14.5 months (SD 27.3), 34 patients had a dead status (47%). Based on response to chemotherapy drug, 41 of patients have good prognosis which overall survival beyond 36 months, other 31 patients have poor prognosis, which had been dead in 18 months.

Identification of somatic mutations in Gastric Carcinomas
In all of the 72 primary tumor samples an adequate library for sequencing was obtained. After read quality filtering, mapping, and alignment to the reference genome (hg19), the median sequencing depth was 277×. A total of 17993 genetic variations (median 207, range 14-959, SD 107.7) were identified in exonic (including nonsense, missense, nonstop, frame_shift_del, frame_shift_ins, in_frame_del and in_frame_ins) ( Figure 1A). The predominant nucleotide changes were cytosine to thymine transitions, which is consistent with previous gastric carcinomas genomics data ( Figure 1B).
The mean number of mutations was 150, range 14-959, SD 107.7. These variants were distributed across 1495 different genes and were found in 72 samples (92% with at least two mutations). 4985 of the missense mutations identified are predicted to be possibly or probably damaging by Sift and Polyphen2, as these are more likely to play a causative role in tumorigenesis and represent therapeutic vulnerabilities. The number of damaging mutations considered to be chemo_sensitive samples was 71 (range 7-292), and 67 in chemo_resistant samples (range 6-233). The somatic mutation rates were similar for chemo-sensitive and chemo-resistant (average mutation rate of 243 versus 254 mutations per tumor, p = 0.79). Mutations in TTN, MUC16 and FSIP2 were the most frequently observed variants in this study (93%, 71%, and 69%, respectively). No significant difference between groups in Lauren type, vascular invasion and perineuronal invasion (Supplemental Figure 1).

The mutation burden and the Clonality do not correlate with chemotherapy
It has been reported that tumour mutational load or tumour mutational burden was associated with response to immunotherapy in multiple cancer types. To examine whether it also correlated with the chemotherapy in patients with GC, we investigated the number of missense mutations as tumor mutation burden in this tumor type. The number of missense mutations in these patients ranged from 13 to 794 (median value, 203). There were no significant differences between groups in chemotherapy (p = 0.34). We next investigated whether there was an association between clonality and chemotherapy. Clonality was positively correlated with tumor heterogeneity, as determined by MATH analysis. There had no significant statistical differences between groups in chemotherapy (p = 0.92). We next performed a survival analysis to evaluate the clinical relevance of clonality, and found no significant association between OS and mutation burden or clonality (Supplemental Figure 2). These results revealed that the number of missense mutations and clonality failed to correlate significantly with chemotherapy survival in GC.

Mutational status and overall survival
In this study, 31 genes had associated with chemotherapy. Then Univariate Kaplan-Meier analysis was used to assess associations between mutations in chemotherapy genes and overall survival from diagnosis. We found a significant association between OS and 16 genes (p < 0.05) (Figure 3 and Figure 4). This analysis revealed that the risk of death was significantly lower for patients with mutant type than for those with wild type, including ABCB5, CACNA1B, CACNA1C, GAL3ST1, GREB1, MICAL2, PTRF and ZNF540 (Figure 3). The other eight genes had high risk of death in mutant type than wild type, including BMPR2, CNTRL, DIDO1, KRIT1, MSH2, NKTR, PRKCB and TCHH (Figure 4).

Discussion
Multiple molecular mechanisms have been reported in drug resistance, such as drug inactivation, drug target alteration, drug efflux, DNA damage repair, cell death inhibition and the epithelialmesenchymal transition (EMT) [14]. However, the specific mechanisms mediating escape from chemotherapy and gastric cancer recurrence remain largely unknown. To address this issue, we target-sequenced a cohort of 72 gastric cancer tissues to find chemotherapy predictive biomarkers and new treatment targets. We found that the commonly recurring genomic aberrations involved TP53, ARID1A, RHOA, RNF43, MUC6, KRAS, PTEN, CDH1, SMAD4, GLI3, ZIC4, ERBB2/3/4, NRG1, DCLK1, and PIK3CA had similar mutation rate between the chemo sensitive subgroup and the chemo resistant subgroup. As a result, sixteen genes including ABCB5, CACNA1B, CACNA1C, GAL3ST1, GREB1, MICAL2, PTRF, ZNF540, BMPR2, CNTRL, DIDO1, KRIT1, MSH2, NKTR, PRKCB and TCHH had significantly difference between the two groups ( Figure 5).  Firstly, activation and inactivation of cancer chemotherapeutic agents can play a major role in the development of resistance. One example of this is observed in the treatment of gastric cancer, capecitabine that is activated after phosphorylation that convert it to 5-fluorouracil [15]. However, hypermethylation of the promoter of gene encoding thymidine phosphorylase is causing capecitabine resistance [16]. Secondly, the damaged DNA repair has also a clear role in anticancer drug resistance [17]. DNA mismatch repair (MMR) in the DNA damage response mediates the cytotoxicity of some DNA-damaging anticancer agents; for example, hypermethylation of MLH1 has been reported to cause resistance to cisplatin and carboplatin [18]; and the absence of MSH2 has been shown to increase resistance to thiopurines [19]. But in this study, the mutation of MSH2 was identified to increase sensitive to chemotherapy, with the mutation rate was 24%. Thirdly, Genomic instability that is called chromosomal instability (CIN) drives cancer development and increases chemotherapy resistance [20,21]. The CIN gene (NEK2) mediated drug resistance was further demonstrated to be through AKT-mediated upregulation of ABC transporters [22]. Several studies have clearly implicated centrosome has a role in tumorigenesis and centrosome dysfunction can result in chromosome instability [23]. Centriolin is encoded by CNTRL gene and involves in maintaining centrosome structure and functions [24]. In our study, the higher mutation rates of gene CNTRL in chemo_sensitive subgroup ensure that they have a decreased risk of chromosomal instability in tumorigenesis.
Additionally, PRKCB, CACNA1A/B and BMP signaling pathway mainly function in regulating cell proliferation. The gene PRKCB may activate the MAPK pathway, regulating cell proliferation [25]. And a substantial alteration of Ca2+ could be responsible for cell proliferation and cell apoptosis through the cell cycle [26]. Bone morphogenetic protein (BMPs) is a large family of growth factors that belongs to tumor growth factor-β superfamily. They have crucial roles in regulating a wide range of developmental functions, such as pattern formation, differentiation, proliferation, migration and so on [27]. BMP protein binding to type 1 and 2 receptors (BMPR1 and BMPR2), activated BMP signaling, increased migration and invasion. In gastric cancer, bone morphogenetic protein-2 (BMP2) significantly inhibited cancer cell proliferation by suppressing the expression of β-catenin, c-Myc, and aurora kinases [28]. Death inducer-obliterator 1 (DIDOI) as a BMP-specific Smad-regulated target gene, promotes the attachment, migration, invasion and apoptosis resistance of melanoma cells [29]. In our study, there is a significant difference in the mutation rate of PRKCB, CACNA1A/B, BMPR2, and DIDOI in different groups.
The tumor microenvironment (TME) not only plays a vital role during cancer progression and metastasis but also has profound effects on therapeutic efficacy. In our study, mutation rates among genes related to cytoskeleton (MICAL2) [30], cell adhesion (KRIT1) [31,32], were also significantly difference (MICAL2 p = 0.005, KRIT1 p = 0.018).
Several the ATP-binding cassette (ABC) transporter family of transmembrane proteins such as ABCB1 (MDR1, P-glycoprotein), ABCC1 (MRP1, multidrug resistance-associated protein 1) and ABCG2 (BCRP, Breast Cancer Resistance Protein), has been linked to mediate multidrug resistance of human cancer by promoting anticancer drug efflux [33,34]. And ABCB5 (ATP-binding cassette member B5) was also as a clinically relevant multidrug resistance mediator in colorectal cancer [35], hepatocellular carcinoma [36] and human malignant melanoma [37]. In this study we observed that 23 % of the tumors had an ABCB5 mutation in chemoresistant subgroup and 5 % in chemosensitive subgroup, indicating that the mutation of ABCB5 plays an important role in drug resistant in gastric cancer. And PTRF knockdown also reduced multidrug resistance in cancer cells via the fortification of lipid rafts [38]. In this study, the mutant PTRF increased the drug resistance in gastric cancer.

Conclusion
Two important conclusions can be inferred from these data. First, genes including TP53, CTNNB1, ARID1A, CDH1, PTEN, KRAS and so on, which were as drivers of gastric carcinogenesis, play an important role in tumor development and progression. However most of these genetic variants do not influence the drug responses. Second, drug resistance to chemotherapy in gastric carcinomas is a complex and integrated activation and inactivation, changes of apoptosis and proliferation, drug efflux, DNA damage repair, and the tumor microenvironment. To discover and prioritize recurrent mutations associated with response to chemotherapy agents, larger sample series are urgently needed.