Noam Goldberg
email: noam.goldberg@biu.ac.il

Michael Poss
email: michael.poss@lirmm.fr

Linearized Formulations for Failure Aware Barter Exchange

Keywords: Longest path, MILP, layered formulation, clearing algorithms, kidney exchange, failure aware barter exchange

 this NP-hard problem was solved heuristically or exactly only for limited path lengths. Although the problem appears highly nonlinear, we formulate it as a mixed-integer linear program (MILP). A computationally tractable layered formulation that approximately solves larger instances is also proposed and a computational study is presented for evaluating the proposed formulations.

Introduction

In this paper we consider a maximum expected value partially successful (s, t)path problem (PSPP), where each (s, v)-subpath of an (s, t)-path, up to some vertex v along the path, has a value of the product of probabilities of arcs along the (s, v)-subpath, times the sum of arc profits.

Kidney exchange market clearing problems are barter exchanges where donors and recipients interact without exchange of monetary means; pairs of recipients and donors enter the market when they are not physiologically compatible and a suitable exchange is sought so that the recipient can be accommodated. Paired-kidney donation graphs indicate donor-recipient compatibility using directed arcs. Each vertex in such a graph corresponds to a pair of a recipient (patient) and donor -typically a family member who is willing to donate a kidney in exchange for the related patient receiving one [START_REF] Abraham | Clearing algorithms for barter exchange markets: enabling nationwide kidney exchanges[END_REF][START_REF] Rees | A nonsimultaneous, extended, altruistic-donor chain[END_REF][START_REF] Roth | Efficient kidney exchange: Coincidence of wants in markets with compatibility-based preferences[END_REF]. The market clearing problem is that of matching all recipients with donors via cycles [START_REF] Roth | Efficient kidney exchange: Coincidence of wants in markets with compatibility-based preferences[END_REF], and paths [START_REF] Abraham | Clearing algorithms for barter exchange markets: enabling nationwide kidney exchanges[END_REF]; see also later work [START_REF] Constantino | New insights on integerprogramming models for the kidney exchange problem[END_REF][START_REF] Anderson | Finding long chains in kidney exchange using the traveling salesman problem[END_REF] and references therein. Paths or chains can be formed when altruistic donor vertex is available to initiate the exchange chain without requiring any donation in return (there is no related recipient and so an altruistic vertex by definition only has outgoing arcs in the donation graph). Exchange cycles and chains are illustrated using small compatibility graph examples in Figure 1. Additional work concerning kidney-transplant exchanges using cycles and/or paths initiated by altruistic donors includes [START_REF] Glorie | Kidney exchange with long chains: An efficient pricing algorithm for clearing barter exchanges with branchand-price[END_REF][START_REF] Manlove | Paired and altruistic kidney donation in the UK: algorithms and experimentation[END_REF].

There is significant uncertainty with respect to the success of a particular donation to take place as a part of the exchange: Physiological incompatibilities may be detected at later stages prior to the transplantation. In addition, a particular donor may become unavailable due to illness or simply because reneging on the agreement. Yet, chains are considered less risky than cycles because a failure of one donation does not necessarily prevent the entire chain of transplants from taking place. This setting with uncertainty gives rise to partially successful paths and the objective considered in this paper; see the expected utility of chains formulated in [START_REF] Dickerson | Failure-aware kidney exchange[END_REF] and the future work discussion in [START_REF] Li | Optimal decisions for organ exchanges in a kidney paired donation program[END_REF]. The corresponding problem for a cycle, where all edges must succeed in order for it to be deemed successful, is also considered in [START_REF] Dickerson | Failure-aware kidney exchange[END_REF] as well as other work including [START_REF] Chen | Graphbased optimization algorithm and software on kidney exchanges[END_REF].

Most literature to date has considered heuristic methods for solving the pricing subproblem of barter-exchange market clearing optimization problem. Yet some exact solutions approaches have been considered in the absence of uncertainty, as [START_REF] Anderson | Finding long chains in kidney exchange using the traveling salesman problem[END_REF][START_REF] Glorie | Kidney exchange with long chains: An efficient pricing algorithm for clearing barter exchanges with branchand-price[END_REF] that consider a cycle and path subproblems (i.e., matching failures). On the other hand, the pricing problem of the failure-aware market clearing problem introduced in [START_REF] Dickerson | Failure-aware kidney exchange[END_REF], has been solved by heuristic means for both cycles and paths. In particular, it appears that exact approaches for the failure-aware path subproblem (i.e., PSPP) has not been considered in previous work. However, solving the pricing problem exactly is important for two reasons. First, exact column generation and/or branch-and-price schemes ultimately need to solve the pricing problem exactly at least in the last few iterations. Second, PSPP is arguably of interest in its own right, and more so than the corresponding cycle problem, because as altruistic donors are scarce and their offer to donate may expire, the PSPP may arise in an online setting.

Uncertainty of the planned donations and failed matches are also considered in stochastic setting with recourse where the probabilities of success are known (similar to the probabilities given as a part of the problem's input in the current setting), for example in [START_REF] Klimentova | Maximising expectation of the number of transplants in kidney exchange programmes[END_REF], or in adjustable robust settings with finite uncertainty sets [START_REF] Carvalho | Robust models for the kidney exchange problem[END_REF]. Both stochastic and robust approaches consider different types of recourse such as "back-arcs recourse" [START_REF] Klimentova | Maximising expectation of the number of transplants in kidney exchange programmes[END_REF][START_REF] Carvalho | Robust models for the kidney exchange problem[END_REF], where exchanges Here compatibility is assumed to depend only on the blood types but in practice compatibility graphs may be sparser due to additional compatibility considerations. may take place along the planned chains (or cycles) up to the failed match. This is in contrast with complete recourse that allows for re-planning of entire cycles and chains in the event of failure. The current model with kidney donation chains maximizes the expected value while allowing for such "back-arcs" recourse.

Our problem bears similarities with the probabilistic all-or-nothing problem, which is to select a subset of items so as to maximize the expected total profit -the product of probabilities of selected items times the sum of profits of these items. In particular, the all-or-nothing probabilistic path problem [START_REF] Goldberg | Maximum probabilistic all-or-nothing paths[END_REF], has the subsets of items corresponding to paths in a graph, each of which has an expected value that is equal to the product of its underlying arc probabilities times the sum of its arc profits. Other all-or-nothing problems have been considered over all unrestricted subsets of a given ground set, and matchings in a graph [START_REF] Goldberg | On the complexity and approximation of the maximum expected value all-or-nothing subset[END_REF]. In the context of failure-aware barter-exchange this "all-ornothing" objective is considered more appropriate for cycles than for paths [START_REF] Dickerson | Failure-aware kidney exchange[END_REF].

Next we formally describe PSPP and relate our mixed-integer nonlinear programming (MINLP) formulation to the optimization problem as it appears in the literature. We then propose a linearization and reformulation of the MINLP, followed by an extended formulation based on a layered graph. Finally, we conclude with computational experiments on simulated kidney exchange data to evaluate the proposed formulations along with attempts to strengthen them.

Problem definition

We now formally state the (PSPP) problem. Let G = (V, E) be a directed graph, s, t ∈ V be the source and the sink, and let δ + (u) = {v ∈ V : (u, v) ∈ E} and δ -(u) = {v ∈ V : (v, u) ∈ E} be sets of direct successors and predecessors of vertex i ∈ V , respectively. The problem with multiple sources and sinks can be also handled, by constructing a supersource s and a supersink t that are connected to the multiple sources and sinks, respectively. For each edge e = (u, v) ∈ E, we are given a positive profit c e = c uv ≥ 0 and a probability of success p e = p uv ∈ [0, 1]. Let m = |E| and n = |V |. Given a path π ∈ P(t) from s to t, where P(t) ∈ 2 m is the set of all elementary paths from s to t (meaning every vertex appears at most once), for v ∈ π we define π(v) as the subpath of π joining s and v. Following [START_REF] Dickerson | Failure-aware kidney exchange[END_REF] (see also the future work discussion in [START_REF] Li | Optimal decisions for organ exchanges in a kidney paired donation program[END_REF]), the maximum expected value PSPP is to find a path π ∈ P(t) so that the objective function

z(π) = (u,v)∈π (1 -p uv ) e∈π(u) c e f ∈π(u) p f + e∈π c e f ∈π p f , (1)
is maximized. First, the following proposition simplifies objective (1) to a form that can be exploited in order to propose an MILP for this problem.

Proposition 1 The objective (1) is equivalent to

z(π) = (u,v)∈π c uv f ∈π(v) p f (2)
Proof By straightforward algebra,

p f + (u,v)∈π e∈π(u) c e f ∈π(u) p f + e∈π c e f ∈π p f = - (u,v)∈π e∈π(u) c e f ∈π(v) p f + (u,v)∈π e∈π(v) c e f ∈π(v) p f = (u,v)∈π c uv f ∈π(v) p f = z(π)
Then, the maximum expected value PSPP is max{z(π) :

π ∈ P(t)}. (3)
Evidently the problem (3) with p e = 1 for all e ∈ E is a longest path problem that is strongly NP-hard in general directed graphs [START_REF] Garey | Computers and Intractability: A Guide to the Theory of NP-Completeness[END_REF].

Note that although a particular source and destination are assumed there can be many altruistic donors and a chain may not have to end in any particular recipient. Accordingly, a super-source s is created with arcs (s, a) having

s a 1 v 1 v 2 v 3 t v 4 v 5 v 6 a 2 0, 1 0, 1
probability p sa = 1 and profit c sa = 0 for all altruistic donors a ∈ V . Similarly the super-sink t with arcs (v, t) having probability p vt = 1 and c vt = 0 reflects the fact that any recipient may be last in the transplant chain. An example of such a graph with two altruistic donors and six donor-recipient pairs is shown in Figure 2. Finally, note that although altruistic donors are typically used to initiate chains, universal recipients (i.e., patients with an AB blood type) corresponding to vertices that have an in-degree of n -3, may also be effectively used, since they are compatible with each and every possible donor who may terminate the exchange.

MILP formulations

We now present mixed-integer linear programming (MILP) formulations to exactly and approximately solve (3).

Exact formulation

Although the problem (3) appears to be highly nonlinear in fact it can be formulated as a MILP. Given a graph G = (V, E), with s, t ∈ V , let the path flow polytope be denoted by

Π = x ∈ R |E| + u∈δ -(v) xuv= u∈δ + (v) xvu for v ∈ V \ {s, t} u∈δ + (s) xsu= u∈δ -(t) xut=1 . ( 4
)
and consider the following MILP formulation. max

(u,v)∈E c uv q uv (5a)
subject to:

w∈δ -(u)

q wu = v∈δ + (u) 1 p uv q uv u ∈ V \ {s, t} (5b)
q uv ≤ M uv x uv (u, v) ∈ E (5c) x ∈ Π (u, v) ∈ E (5d) x uv ∈ {0, 1}, q uv ≥ 0 (u, v) ∈ E, (5e)
where 0 < M uv ≤ 1 is a sufficiently large probability constant. For each (u, v) ∈ E the decision variable x uv is used to indicate whether (u, v) is selected as a path edge, and the auxiliary variables q uv is used to define q uv = (u,v)∈π(v) p uv as in z(π). The constraints corresponding to (5d) are the s -t path-flow conservation constraints. Constraint (5c) enforces that q uv = 0 for each (u, v) ∈ E with x uv = 0. Finally, (5b) implements the recursion to accumulate the probability product for each pair of selected incident edges (w, u) and (u, v) with x uv = x wu = 1, it requires that q wu p uv = q uv . To strengthen the formulation (5) it is desirable to select M uv for each (u, v) ∈ E as small as possible. It can be observed it suffices to set M uv = p(u)p uv where p(v) denotes the maximum probability of a path from s to v for all v ∈ V \ {s} and p(s) = 1. It is known that p can be determined in polynomial time by solving a maximum reliability problem. In particular, it can be observed by considering the shortest path problem with arc costs -ln p e for each e ∈ E (in an exact arithmetic setting) or replacing the addition in the Bellman equations with multiplication.

To prevent solutions that are feasible for (5) from containing cycles, cyclebreaking inequalities must be added. We consider two alternative ways of breaking the cycles. First, following [START_REF] Taccari | Integer programming formulations for the elementary shortest path problem[END_REF], we consider generalized cutset inequalities,

(u,v)∈δ + (S)
x uv ≥ (w,v)∈δ + (w)

x wv S ⊆ V \ {s, t} : |S| ≥ 2, w ∈ S. (6)
These constraints break any cycle by enforcing that the number of arcs leaving S be not smaller than the number of arcs outgoing from any vertex k of S.

A second type of inequalities will prevent cycles if the probabilities are strictly less than unity (specifically their product along cycles). The sufficiency of the inequalities

v∈δ + (u) x uv ≤ 1 u ∈ V \ {t}, (7)
for prevention of cycles, given the condition that success is uncertain along cycles, is established by the following proposition.

Proposition 2

The problem (5) with (7) has an acyclic optimal solution if

f ∈c p f < 1 for each cycle c ⊆ E.
Proof Consider (x, q) that is a feasible solution of ( 5) with [START_REF] Constantino | New insights on integerprogramming models for the kidney exchange problem[END_REF]. Note that by the flow-conservation constraints and the integrality of x, the support of x is a union of one path π from s to t and a set of cycles C ⊂ 2 E . Suppose that there is some u ∈ V that is incident to edges in two different cycles, or an edge in a cycle and an edge in π. Then, by the integrality of x, v∈δ + (u) x uv > 1, thereby violating [START_REF] Constantino | New insights on integerprogramming models for the kidney exchange problem[END_REF]. It follows that the cycles c ∈ C and the path π are all vertex disjoint. Now further suppose for the sake of deriving a contradiction that (x, q) is an optimal solution with minimal |C| > 0 and consider a cycle c ∈ C. If q e = 0 for each e ∈ c, then a new solution (x , q) can be defined (using (x, q)) by setting x e = 0 for each e ∈ c. The objective value of (x , q) is equal to that of (x, q) and has one less cycle, thereby establishing a contradiction. Otherwise, if q e > 0 for some e ∈ c, constraints (5b) written for all (u, v) ∈ c imply that f ∈c p f = 1.

Rounded formulation

The formulation from the previous section tends to have a weak linear programming relaxation, due the big-M coefficients involved in constraints (5c). We present in this section an alternative and stronger approach that is based on layered graphs, where multiple copies of graph vertices and edges are created and "stacked on top of each other"; see for example [START_REF] Gouveia | Layered graph approaches for combinatorial optimization problems[END_REF]. Specifically, we consider next a base 0 < α < 1 and define the rounded logarithms of the probabilities as r e = log α p e . The set of rounded logarithms values is K = {0, . . . , n max e∈E {r e }}. Define V 0 = V \ {s, t} and let its induced subset of edges be denoted by E 0 = E[V 0 ]. Then, consider copies of V 0 , one for each k ∈ K, to be denoted by

V k 0 = {(v, q) ∈ V 0 × K | q = k }; each mem- ber (v, k) ∈ V k 0 is the k-th copy of v ∈ V 0 in V . The layered graph is then G = ( V , E), where V = {(s, 0), (t, 0)} ∪ k∈K V k 0 and E = (u, k, v, k + r u,v ) (u, v) ∈ E 0 , k ∈ {0, . . . , n max e∈E {r e } -r u,v } ∪ (s, 0, v, r s,v ) v ∈ δ + (s) ∪ (u, k, t, 0) u ∈ δ -(t), k ∈ K .
Notice that G is acyclic. Let Π be the corresponding path flow polytope Π given by ( 4), defined for the layered graph G in place of G, with (s, 0) and (t, 0) as the source and the sink vertices. Evidently, all paths in this graph correspond to extreme points of this polytope. Consider the formulation max

e=(u,k,v, k)∈ E α kc uv x e (8a)
subject to:

v∈V \{s}:∃ k∈K e=(u,k,v, k)∈ E x e ≤ 1 u ∈ V \ {t} (8b)
x ∈ Π (8c)

x e ∈ {0, 1} e ∈ E. (8d)
Note that while any path π = ((u 1 , q 1 ) = (s 1 , q 1 ), (u 2 , q 2 ), . . . , (u l , q l ) = (t l , q l )) in G is elementary, its projection onto V , π = (u 1 = s, u 2 , . . . , u l = t) might contain a cycle. We forbid such paths π (whose projection contain one or more cycles) by limiting the number of outgoing arcs from any u ∈ V \ {s} to at most one in constraints (8b). Also note that the rounding makes solving formulation [START_REF] Dickerson | Optimizing kidney exchange with transplant chains: Theory and reality[END_REF] an inexact solution of (5) (and ( 3)). However, if α is chosen close enough to one, the worst-case performance of the optimal solution of ( 8) can be bounded from the theoretical viewpoint as formalized below. 1) , then each solution that is optimal to (8) is an (1 -)-approximate solution to (5) with [START_REF] Constantino | New insights on integerprogramming models for the kidney exchange problem[END_REF].

Observation 1 For α ≤ (1 -) 1/(n-
The correctness of this observation follows from the fact that if path π = (u 1 , . . . , u l ) is optimal in the original graph then there is a path π in G, whose projection onto V is π, and whose objective value is

(u,k,v, k)∈ π α kc uv = l-1 i=1 α i j=1 ru j u j+1 c uiui+1 ≥ l-1 i=1 c uiui+1 α i i j=1 p uj uj+1 ≥ α n-1 z(π).

Computational results

In this section we experiment with simulated kidney compatibility graphs that are publicly available [START_REF] Mattei | Preflib: A library for preferences http://www. preflib. org[END_REF] or randomly generated following the structure suggested in kidney exchange literature. All data used is maintained in a public repository [14].

Graph topology simulated according to [START_REF] Ashlagi | Nonsimultaneous chains and dominos in kidney-paired donation-revisited[END_REF] (graph[n]a[A] datasets). We generate a "sparse-dense" graph with the vertex connectivity generated according to the panel reactive antibodies (PRA) levels as suggested in [START_REF] Ashlagi | Nonsimultaneous chains and dominos in kidney-paired donation-revisited[END_REF], and we make use of additional data given in [START_REF] Ashlagi | The need for (long) chains in kidney exchange[END_REF][START_REF] Dickerson | Failure-aware kidney exchange[END_REF]. In particular, graph vertices are created with a highly sensitized recipient with probability 0.27. Such highly sensitized patients (with a high PRA value) are less likely to be compatible with a random donor. Accordingly, these recipients are then connected with donors (corresponding to an incident incoming edge) with probability 0.03. The failure probabilities are assigned also based on PRA levels and highly sensitized patients are matched with a success probability of 0.5. Otherwise, low PRA recipients are compatible with donors with probability 0.5, 57% of them are matched with a success probability of 0.75 (cross-match failure rate of 0.25) and 43% of them are matched with a success probability of 0.95 (cross-match failure rate of 0.05). The values of matching PRA recipients are generated from a normal distribution with mean 10 and standard deviation 2. The values of matching highly sensitized recipients are a factor of 1.5 greater (so the mean is 15) as suggested in [START_REF] Dickerson | Failure-aware kidney exchange[END_REF].

Graphs based on simulated donor data instances of [START_REF] Dickerson | Optimizing kidney exchange with transplant chains: Theory and reality[END_REF] and [START_REF] Dickerson | Failure-aware kidney exchange[END_REF] (MD[n] datasets). These instances were generated using a simulator based on [START_REF] Saidman | Increasing the opportunity of live kidney donation by matching for two-and three-way exchanges[END_REF]. Graph arcs were constructed based on blood-type compatibility and we set the weights to be either unitary in datasets MD[n]unit or normally distributed with mean 10 and standard deviation 2 in datasets MD[n]-stoch. For patients that are considered highly sensitized with PRA levels exceeding 0.73 in the original data, the arc values considered were twice as much as initially generated. The success probabilities are set to 0.5 for highly sensitized patients (determined by the PRA levels). Otherwise, a success probability of 0.75 is set with probability 0.57 and success probability of 0.95 is set with probability 0.43 (similar to the graph[n]a[A] datasets).

Implementation details

All algorithms have been coded in Julia 1.3.1 on a platform using a (Intel i7-10510U) 1.80GHz CPU with 4 cores and an 8MB cache and 16GB of RAM. The mathematical programs are modeled using the JuMP v0.21.1 [START_REF] Dunning | JuMP: A modeling language for mathematical optimization[END_REF] and Gurobi v0.7.6 packages and solved with the Gurobi 9.0.1 solver running 4 threads. None of our instances contain an edge with a probability equal to 1. Although, if an instance has multiple altruistic donors, we then connect supersources and supersinks to all vertices with edges having probability 1, but this may not lead to cycles c such that f ∈c p f = 1. Hence, the condition of Proposition 2 holds so inequalities [START_REF] Constantino | New insights on integerprogramming models for the kidney exchange problem[END_REF] are enough to guarantee the existence of optimal acyclic solutions. This being said, inequalities [START_REF] Chen | Graphbased optimization algorithm and software on kidney exchanges[END_REF] are not dominated by inequalities [START_REF] Constantino | New insights on integerprogramming models for the kidney exchange problem[END_REF], so separation of the former through callbacks might further improve the overall performance of the formulation.

Following [START_REF] Taccari | Integer programming formulations for the elementary shortest path problem[END_REF], we thus tested the separation of inequalities (6) at both fractional and integer solutions by checking whether they are satisfied on each strongly connected component of the graph induced by the positive components of x that is optimal to (a relaxation of) [START_REF] Carvalho | Robust models for the kidney exchange problem[END_REF]. Our results indicated a worsening of the overall performance: the relaxation improvement, often nonexistent, is too little to compensate for the time spent in the separation procedure. Therefore, we do not use inequalities [START_REF] Chen | Graphbased optimization algorithm and software on kidney exchanges[END_REF] in the results presented next.

Description of the results

Table 1 shows for the particular kidney exchange graph instances that we experimented, statistics and characteristics including the number of vertices, edges, number of altruistic donors, length of the longest path from the supersource to any vertices in the graph, and the root relaxation gap (in %). Solution times are reported in seconds in Tables 2 and3. Here T indicates that a time limit of 1800 seconds is reached and M indicates that the solver has exhausted the available memory. Several observations can be made from the results of these tables. First, we see that the extended formulation [START_REF] Dickerson | Optimizing kidney exchange with transplant chains: Theory and reality[END_REF] scales reasonably well, as long as the memory limit is not exceeded. Some of the large MD instances, namely 84unit, 84stoch, 125unit, are solved faster than formulation [START_REF] Carvalho | Robust models for the kidney exchange problem[END_REF]. Second, the exact formulation [START_REF] Carvalho | Robust models for the kidney exchange problem[END_REF] appears to be more efficiently solved with instances MD than with instances "graph". To further investigate this behavior we refer to the graph instance characteristics in Table 1. The table data shows that graph30a2 and instances MD43 -MD70 have comparable numbers of vertices and edges. However, graph30a2 cannot be solved in 1800 seconds, while these MD instances are all solved within a few seconds, except instance 66unit which requires roughly 72 seconds. A partial explanation of this behavior could lie in the strong root relaxation gap of ( 5) for these MD instances, ranging from 6 to 20% while for the graph30a2 instance it equals 31%. The strength of the relaxation appears to correlate with the length of paths in these graphs (directly impacting M uv ); for the MD instances with 34 vertices the length ranges from 17 to 22 hops, compared with 30 hops in graph30a2.

Next, the third group of gap columns of Tables 2 and3 provide the ratio of the best solution found by the extended formulation [START_REF] Dickerson | Optimizing kidney exchange with transplant chains: Theory and reality[END_REF] for different values of α, to the best solution determined by the exact formulation (5) (expressed as %). All instances in this table may exceed the time limit when solved exactly by formulation [START_REF] Carvalho | Robust models for the kidney exchange problem[END_REF] and thus lead to negative gap values. Evidently, α < 0.8 seems to lead to poor solutions, while setting α higher than 0.8 appears to yield only marginal improvements. Finally, we remark that for the instance MD125unit that cannot be solved to optimality, the layered formulation (8) even provides slightly better solutions than the best solution found by the exact formulation (within the 30 minute time limit). Note that the theoretical gap of Observation 1 for entries in these tables ranges from ≤ 1-α n-1 ≤ 43% for n = 12, α = 0.95 to more than 99% for most table entries where α ≤ 0.9 and n ≥ 50 (or even n ≥ 20 when α = 0.75). Evidently, the gap appearing in practice tends to be much smaller than the theoretical guarantee.

Conclusions and future work

In this paper we proposed a linear formulation for the maximum expected value partially successful path problem. For this formulation we prove and evaluate cycle prevention inequalities that outperform standard longest path cycle prevention inequalities. We then evaluated a layered formulation that effectively solved the problem within a guaranteed approximation gap for larger instances. According to our experiments the approximation gap in practice usually appeared much smaller than the theoretical guarantee. The problem considered in the current paper is strongly NP-hard in general graphs. In contrast, computing a path that maximizes the similar yet different all-or-nothing objective remains NP-hard even in directed acyclic graphs

  Example of a chain initiated by an altruistic donor.

Fig. 1

 1 Fig.1Examples of a cycle and chain in organ compatibility graphs. Blood types of recipients and donors are indicated on the top part and lower part, respectively, of each vertex. Dashed and solid edges indicate unmatched or matched donations, respectively. Here compatibility is assumed to depend only on the blood types but in practice compatibility graphs may be sparser due to additional compatibility considerations.

Table 1

 1 Kidney exchange graph instance main characteristics. LP stands for longest path and the root gaps are in % with respect to formulation[START_REF] Carvalho | Robust models for the kidney exchange problem[END_REF]. Here na denotes the number of altruistic donors.

	Instance	|V | |E| na LP Root gap
	MD12stoch	18	96 1	9		9
	MD12unit	18	96 1	9		10
	MD19stoch	18 150 1 14		7
	MD19unit	18 150 1 14		4
	MD21stoch	18 114 2 10		7
	MD21unit	18 114 2 10		1
	MD23stoch	18 121 2 13		17
	MD23unit	18 121 2 13		28
	MD43stoch	34 402 1 19		14
	MD43unit MD44stoch MD44unit MD51stoch MD51unit MD60stoch MD60unit MD61stoch MD61unit	34 402 1 19 34 337 1 17 34 337 1 17 34 405 3 19 34 405 3 19 34 371 3 19 34 371 3 19 34 422 4 19 34 422 4 19		13 15 11 6 16 10 9 16 17	Instance graph20a1 22 178 1 20 |V | |E| na LP Root gap 25 graph30a2 32 329 2 30 31 graph40a2 42 681 2 41 30 graph50a3 52 981 3 51 23 graph60a3 62 1387 3 61 17 graph70a4 72 1851 4 71 13 graph80a4 82 2247 4 81 10
	MD66stoch	34 403 4 22		20
	MD66unit	34 403 4 22		16
	MD70stoch	34 332 4 18		13
	MD70unit	34 332 4 18		18
	MD81stoch	66 2205 3 55		15
	MD81unit	66 2205 3 55		20
	MD84stoch	66 1512 3 36		29
	MD84unit	66 1512 3 36		27
	MD125stoch 130 6026 6 73		7
	MD125unit 130 6026 6 73		12
				CPU Seconds				Gap (%)
	Instance		α 0.75	0.8 0.85	0.9 0.95 (5) 0.75 0.8 0.85 0.9 0.95
	graph20a1		3	6	7	10	16 11	31	4	4	1	0
	graph30a2		14	32	40	79	25	T	44	3	1	0	0
	graph40a2		119	345	635	115 608	T	55	1	0 -1	-1
	graph50a3		1198	128	407	590 114	T	49	0	-1 -1	-1
	graph60a3		233	386	534 1456 216	T	71	0	-1 -1	-1
	graph70a4		T 1521 1784	729 662	T	72	0	0 -1	-1
	graph80a4		T	T	T	492	T	T	74	1	20	0	-

Table 2

 2 Solution CPU seconds and and optimality gap for approximate MILP (8) for the "graph" instances.

ca 1 v 1 , pa 1 v 1 ca 1 v 3 , pa 1 v 3 cv 1 v 2 , pv 1 v 2 cv 2 v 3 , pv 2 v 3 ca 1 v 1 , pa 1 v 1 ca 1 v 1 , pa 1 v 1 ca 1 v 2 , pa 1 v 2 cv 3 v 4 , pv 3 v 4 cv 6 v 4 , pv 6 v 4 cv 3 v 5 , pv 3 v 5 0, 1 0, 1 0, 1 0, 1 cv 5 v 6 , pv 5 v 6 cv 6 v 5 , pv 6 v 5 0, 1 0, 1Fig.2Example organ compatibility graph with two altruistic donors (a 1 and a 2 ) and six vertices corresponding to donor-recipient pairs (v 1 , . . . , v 6 ).

CPU Seconds

Gap (%) Instance α 0.75 0.8 0.85 0.9 0.95 (5) 0.75 0.8 0.85 0.9 0.95

Table 3 Solution CPU seconds and optimality gap for approximate MILP [START_REF] Dickerson | Optimizing kidney exchange with transplant chains: Theory and reality[END_REF] for the MD instances.

(DAGs); see [START_REF] Goldberg | Maximum probabilistic all-or-nothing paths[END_REF]. Although it is not difficult to extend the FPTAS suggested for the all-or-nothing objectives to the current problem, it is a question of theoretical interest for future work to determine whether the current problem remains NP-hard in a DAG. While robust optimization models have been considered where the existence of edges and the value of exchanges are subject to uncertainty [START_REF] Carvalho | Robust models for the kidney exchange problem[END_REF][START_REF] Mcelfresh | Scalable robust kidney exchange[END_REF], the failure-aware model considered in this paper could be extended to a distributionally robust model in which the edge success probabilities are subject to uncertainty. This is also a subject for future work. Finally, computational techniques such as the network layering that is applied to our MIP could potentially be applied to other kidney exchange formulations such as ones appearing in [START_REF] Constantino | New insights on integerprogramming models for the kidney exchange problem[END_REF].


A Missoffe
email: alexia.missoffe@supelec.fr

J Juillard
email: jerome.juillard@supelec.fr.

D Aubry

A reduced-order piecewise-linear model of squeeze-film damping for deformable structures including large displacement effects

Keywords: squeeze-film damping, reduced-order modelling, microswitch, piecewise linear

Introduction

Correct modelling of damping is essential to capture the dynamic behaviour of a MEMS device. Our interest is squeeze-film damping which models the behaviour of a fluid in small gaps between a fixed surface and a structure moving perpendicular to this surface . The lateral dimensions of the surfaces are large compared to the gap and the system is considered isothermal. Squeeze film damping is then governed by the Reynolds equation [START_REF] Bao | Micro mechanical transducers: pressure sensors, accelerometers and gyroscopes[END_REF]:

3 ( ) 12
G GP P P t m ae ö ¶ Ñ Ñ = ç ÷ ¶ è ø , (1)
where ( )

t y x G , ,
is the distance between the moving and the fixed surface, ( )

t y x P , ,
is the pressure, the effective viscosity of the fluid [START_REF] Bao | Micro mechanical transducers: pressure sensors, accelerometers and gyroscopes[END_REF]. For small excitation frequencies or amplitudes the squeezed film behaves as a linear damper. For larger amplitudes or frequencies, the gas has no time to flow away and the pressure builds up creating a stiffening effect coupled to a nonlinear damping. The boundary conditions for (1) are usually chosen as trivial: zero pressure variation or zero pressure gradient , although some authors have considered less ideal, frequency-dependent and aspect-ratio dependent boundary conditions [START_REF] Veijola | Extending the validity of squeezed-film damper models with elongations of surface dimensions[END_REF][START_REF] Pandey | Influence of Boundary Conditions on the Dynamic Characteristics of Squeeze Films in MEMS devices[END_REF]. A complete review on this equation and its different regimes can be found in [START_REF] Bao | Squeeze film air damping in MEMS[END_REF]. Coupling the Reynolds equation to the equation governing the mechanical behaviour of the microstructure leads to a nonlinear system of partial differential equations (PDEs), which has no analytical solution and must be simplified thanks to some assumptions. The most commonly made assumptions are the following:

-uniform displacements, i.e.

0 = ¶ ¶ = ¶ ¶ y G x G
(for example, [START_REF] Veijola | Equivalent-circuit model of the squeezed gas film in a silicon accelerometer[END_REF]).

-steady-state sinusoidal excitation, i.e.

( )

t G G e w sin = [6-7]. -small displacements, i.e. g G G + = 0 and 0 G g <<
, where 0 G is the nominal gap of the structure at rest or close to a static equilibrium [START_REF] Li | A new model for squeeze-film damping of electrically actuated microbeams under the effect of a static deflection[END_REF][START_REF] Younis | Microplate Modeling under Coupled Structural-Fluidic-Electrostatic Forces[END_REF]. -small pressure variations, i.e. P is the ambient pressure. These hypotheses prove to be useful in a variety of applications, if only for gaining insight of nonlinear damping phenomena. However, in many cases, it is difficult to justify their use: for example, it is clear to see that none of the first three hypotheses holds when trying to estimate the switching time of a micro-switch. Most micro-switches do not undergo uniform displacements, nor can these displacements be considered small, and the behaviour of a micro-switch is fundamentally transient.

To date, the most notable attempts to tackle the problem of reduced-order modelling (ROM) of squeeze-film damping with large, non-uniform displacements have been made by Younis et al. [START_REF] Younis | Modeling and Simulation of Microelectromechanical Systems in Multi-Physics Fields[END_REF][START_REF] Younis | Simulation of squeeze-film damping of microplates actuated by large electrostatic load[END_REF][START_REF] Younis | A new approach to the modeling and simulation of flexible microstructures under the effect of squeeze-film damping[END_REF], Mehner et al. [START_REF] Mehner | Reduced-order modeling of fluid-structural interactions in MEMS based on modal projection techniques[END_REF][START_REF] Mehner | Simulation of gas damping in microstructures with nontrivial geometries[END_REF], Yang et al. [START_REF] Yang | Squeeze-film damping for MEMS structures[END_REF][START_REF] Yang | Macromodeling of coupled-domain MEMS devices with electrostatic and electrothermal effects[END_REF], Hung and Senturia [START_REF] Hung | Generating efficient dynamical models for microelectromechanical systems from a few finite-element simulation runs[END_REF] and Rewienski and White [START_REF] Rewienski | A Trajectory Piecewise-Linear Approach to Model Order Reduction of Nonlinear Dynamical Systems[END_REF][START_REF] Rewienski | A trajectory piecewise-linear approach to model order reduction and fast simulation of nonlinear circuits and micromachined devices[END_REF]. In [START_REF] Younis | Modeling and Simulation of Microelectromechanical Systems in Multi-Physics Fields[END_REF][START_REF] Younis | Simulation of squeeze-film damping of microplates actuated by large electrostatic load[END_REF][START_REF] Younis | A new approach to the modeling and simulation of flexible microstructures under the effect of squeeze-film damping[END_REF], the authors propose to solve the nonlinear Euler-Bernoulli beam equation to determine the static deformation of a microplate under a voltage bias. The von Karman plate equations and the compressible Reynolds equation are then linearized close to this operating point and a perturbation method is used to calculate the pressure deviation. In [START_REF] Mehner | Reduced-order modeling of fluid-structural interactions in MEMS based on modal projection techniques[END_REF][START_REF] Mehner | Simulation of gas damping in microstructures with nontrivial geometries[END_REF], the authors use a modal projection method to calculate modal frequency-dependent damping and stiffening coefficients close to a determined operating point. To extend this approach to large displacements, Mehner [START_REF] Mehner | Reduced-order modeling of fluid-structural interactions in MEMS based on modal projection techniques[END_REF] gives an analytical expression of these coefficients as a function of mechanical modal coordinates established by fitting of simulation data for different initial deformations. These approaches are all based on several steady-state sinusoidal calculations [START_REF] Younis | Modeling and Simulation of Microelectromechanical Systems in Multi-Physics Fields[END_REF][START_REF] Younis | Simulation of squeeze-film damping of microplates actuated by large electrostatic load[END_REF][START_REF] Younis | A new approach to the modeling and simulation of flexible microstructures under the effect of squeeze-film damping[END_REF] or simulations [START_REF] Mehner | Reduced-order modeling of fluid-structural interactions in MEMS based on modal projection techniques[END_REF][START_REF] Mehner | Simulation of gas damping in microstructures with nontrivial geometries[END_REF], which increase the for setting up the reduced-order model. The most general approaches may well be those developed in [START_REF] Yang | Squeeze-film damping for MEMS structures[END_REF][START_REF] Yang | Macromodeling of coupled-domain MEMS devices with electrostatic and electrothermal effects[END_REF][START_REF] Hung | Generating efficient dynamical models for microelectromechanical systems from a few finite-element simulation runs[END_REF][START_REF] Rewienski | A Trajectory Piecewise-Linear Approach to Model Order Reduction of Nonlinear Dynamical Systems[END_REF][START_REF] Rewienski | A trajectory piecewise-linear approach to model order reduction and fast simulation of nonlinear circuits and micromachined devices[END_REF]: the authors rely on fully-coupled, nonlinear transient simulations of the complete system (usually a micro-switch) to establish a reduced-order model of the microstructure. These approaches are very general and they can even be successfully applied to the fully nonlinear Reynolds equation [START_REF] Bao | Micro mechanical transducers: pressure sensors, accelerometers and gyroscopes[END_REF]. However, they have a high computational cost (because of the nonlinear/multiphysics/transient simulation they require) and their accuracy depends, to some degree, on the choice of the training trajectory. We present a reduced-order model of the Reynolds equation which is valid for large, nonuniform displacements and transient excitation but with the restriction of small pressure variations. This hypothesis is valid for a microswitch as shown by the results presented in [START_REF] Juillard | A reduced-order model of squeeze-film damping for deformable structures including large displacement effects[END_REF].

Construction of the reduced-order model

We work on the variable 0

P P p - =
, supposing 0 P p << . (1) has then the following form:

÷ ÷ ø ö ç ç è ae ÷ ÷ ø ö ç ç è ae + ¶ ¶ = ÷ ÷ ø ö ç ç è ae Ñ Ñ 0 3 1 12 P p G t p G m (2)
The first reduction step is based on modal projection of (2) which is first transformed via a change of variable on p. The aim of this change of variable is to obtain a spatial operator for which the Laplacian eigenmodes are more relevant than for the operator in (2), conserving its self-adjoint property thus guarantying convergence of the solution. The reduced-order model may be written as:

s x H x f s x A ) ( )) ( ) ( ( = - dt d , (3)
where x and s are respectively the vectors of the mechanical modal coordinates of the moving structure, and of the modes corresponding to the Reynolds equation. For a structure under electrostatic actuation, one may write the full coupled model as:

) ( )) ( ( )) ( ( )) ( ( ( t u t t t dt d z B z f z g + = ( 4
)
where z is the state vector including all modal coordinates. The cost of evaluating the terms g(z(t)), f(z(t)), and B(z(t)) is reduced using a piecewise linear approach described in [START_REF] Rewienski | A Trajectory Piecewise-Linear Approach to Model Order Reduction of Nonlinear Dynamical Systems[END_REF]. The fact that the coefficients in (3) only depend on the mechanical modal coordinates reduces the cost of construction of the piecewise linear model, which has the following structure:

) ( )) ( ) ( )( ( )) ) ( ( ) ( 1 0 1 0 z F z z JF z f z z (z JG z g z i i i elec i s i i i i s i i w w dt d + - + = ú û ù ê ë é - + å å - = - = (5)
where i JG and i JF are respectively the jacobians of the functions g(z) and f(z) at the linearization points i z and (z) F elec corresponds to the electrostatic force. Two problems arise from the piecewise linearization: the choice of the linearization points and the weighting procedure. We choose the linearization points from a simulation of the first reduced model. A new linearization point is chosen when a point is far enough from the already chosen points. We normalize the state variables to calculate the distances so as to take into account all state variables with the most relevance. The weighting procedure is the one described in [START_REF] Rewienski | A Trajectory Piecewise-Linear Approach to Model Order Reduction of Nonlinear Dynamical Systems[END_REF].

Validation of the reduced-order model

We work on the example of a microswitch also treated in [START_REF] Hung | Generating efficient dynamical models for microelectromechanical systems from a few finite-element simulation runs[END_REF][START_REF] Rewienski | A Trajectory Piecewise-Linear Approach to Model Order Reduction of Nonlinear Dynamical Systems[END_REF][START_REF] Rewienski | A trajectory piecewise-linear approach to model order reduction and fast simulation of nonlinear circuits and micromachined devices[END_REF]. We use the piecewise linear reduced order model to determine the switching time of the device for a step voltage between 9 and 10.5 V at atmospheric ambient pressure. The 21 linearization points are chosen along a 9.5V input training trajectory. Fig. 1 shows the response to a 10V input using a linear model, the modal projection model, and the piecewise linearized model. Fig. 2 shows the experimental and simulated data presented in [START_REF] Hung | Generating efficient dynamical models for microelectromechanical systems from a few finite-element simulation runs[END_REF] and results of our piecewise linear reduced model of order 6 for the switching time. Comparison of the experimental and simulated results presented in [START_REF] Hung | Generating efficient dynamical models for microelectromechanical systems from a few finite-element simulation runs[END_REF] to the simulated results obtained with our reduced-order model. The chosen squeeze modes correspond to 1 k =0, 2 and 2 k =1, 3 in (38). For the mechanical part, the first eigenmode is used.

Further work

We can notice that the nonlinear terms depend only on the mechanical modal coordinate and that the model is linear in the modes concerning the Reynolds equation. It can be written as:

z x B x FP z x G ) ( )) ( ) ( ( = - dt d ( 6
)
The piecewise linearization described above doesn t take advantage of this structure. Another piecewise linearized model can be based on the linearization of the terms G, FP and B. As they only depend on the on mechanical coordinate, there is no need for a training trajectory. It is sufficient to discretize in an appropriate way the one dimensional space corresponding to the mechanical coordinate. This is a great advantage as the resulting model does not depend on the relevance of a training trajectory, one the main drawbacks of the model presented above. On the other hand the resulting equation is nonlinear which increases the resolution cost. Full results of this model will be presented.

Figure 1 -

 1 Figure 1-Middle point displacement of a microswitch for a linear model, the projected model and the projected and piecewise linearized model (1 mechanical mode, 4 squeeze modes) for a 10V step voltage. The 21 linearization points are chosen along a 9.5V input trajectory.

Figure 2 -

 2 Figure2-Pull-in time (s) versus applied voltage (V) for 0 P =1.013×105 Pa. Comparison of the experimental and simulated results presented in[START_REF] Hung | Generating efficient dynamical models for microelectromechanical systems from a few finite-element simulation runs[END_REF] to the simulated results obtained with our reduced-order model. The chosen squeeze modes correspond to 1 k =0, 2 and 2 k =1, 3 in (38). For the mechanical part, the first eigenmode is used.


P Soudier

; Faure

; Kushwaha

; Faulon

Cell-Free Biosensors and AI integration

Cell-free biosensors hold a great potential as alternatives for traditional analytical chemistry methods providing low-cost low resource measurement of specific chemicals. However, their large scale use is limited by the complexity of their development.

In this chapter we present a standard methodology based on Computer aided design (CAD) tools that enables fast development of new cell-free biosensors based on target molecule information transduction and reporting through metabolic and genetic layers, respectively.

Such systems can then be repurposed to represent complex computational problems, allowing defined multiplex sensing of various inputs and integration of artificial intelligence in synthetic biological systems.

Introduction :

Metabolite biosensors emerged as a major application of synthetic biology. Repurposing biological systems naturally present in various organisms enabled the development of synthetic sensing devices [START_REF] Van Der Meer | Where microbiology meets microengineering: design and applications of reporter bacteria[END_REF] giving rapid and inexpensive point of care measurement of various molecules of interest levels in a broad range of samples with applications in health [2], environment [3], industries and fundamental research [4]. The potential of these types of devices lies in the fact that they are able to deliver information on the presence and levels of certain chemicals without the need for expensive reagents, trained operators and large size facilities. Transcription factor based biosensors are a subset of those sensing devices repurposing transcription factor and inducible promoter for the quick detection of molecules of interest.

Cell-free biosensors emerged as an alternative for traditional whole-cell biosensors, solving several issues carried by these systems [5]. Cell-free biosensors are efficient, generating big signal over noise ratios [START_REF] Voyvodic | Plug-and-play metabolic transducers expand the chemical detection space of cell-free biosensors[END_REF]. They are suitable for the detection of molecules toxic for bacterial growth or molecules that are not able to cross the cell membrane. Their ability to be freeze dried on paper enables long term storage at room temperature without loss of activity pushing their use as point-of-care devices, especially in low resource communities [6]. The absence of living organisms in it is also facilitating their industrial developments with reduced regulatory issues and low biosafety concerns. Finally the ability to express each gene from distinct DNA fragments (plasmids or linear fragments) [7] and to fine tune their expression levels by varying their respective concentrations enable a fast development and optimisation of new biosensor candidates or more complex devices relying on them. This central property of straightforward tuning of gene expression is necessary for the development of the complex information processing systems that will be described in this chapter.

The main issue encountered in the development of transcription factor (TF) based biosensors is the limited number of molecules for which an interacting transcription factor has been described. The list of these metabolites known to trigger transcriptional response, either in natural systems or in the context of synthetic biosensors, have been described in various databases [8] [9] and compiled in a dataset of small molecules triggering transcriptional and translational cellular responses [10].

In order to detect molecules for which associated transcription factors are not available, new methodologies have been developed relying on the use of metabolic pathways. It has indeed been shown that enzymes could be used to convert non-detectable molecules into molecules known to regulate a characterized transcription factor [11]. This framework has been formalized and showcased as the Sensing-Enabling Metabolic Pathway (SEMP) concept [START_REF] Libis | Expanding Biosensing Abilities through Computer-Aided Design of Metabolic Pathways[END_REF]. This shifted the problem from finding a TF binding the molecule of interest to finding any potential enzymatic route between the molecule of interest and the known set of detectable molecules. This was addressed by repurposing bioinformatic tools traditionally used for retrosynthesis [START_REF] Delépine | RetroPath2.0: A retrosynthesis workflow for metabolic engineers[END_REF] into a platform able to find potential SEMP for any molecule of interest. This web-service called SENSIPATH is a CAD program predicting pathways of up to 2 steps allowing detection for any query compounds [START_REF] Delépine | SensiPath: computer-aided design of sensing-enabling metabolic pathways[END_REF]. This approach was then used for the development of a plug and play cell-free workflow where enzymes formalized as metabolic transducers were used as interchangeable modules converting various chemicals into detectable molecules called effectors [START_REF] Voyvodic | Plug-and-play metabolic transducers expand the chemical detection space of cell-free biosensors[END_REF]. Such biosensors were found to be highly efficient, the cell-free system potentialities allowing to reach optimal response for each biosensor by a systematic tuning of each component involved (reporter, transcription factor and enzyme encoding DNA). They were also applied for the commercially relevant sensing of molecules of interest in complex real-world samples (preservative in beverage, disease biomarker and drug in clinical samples).

In addition to its biosensing potential, this framework has also then been applied to the field of information processing. Indeed, the constructed devices sensing various molecules can also be conceptualised as information processing systems, converting an input signal (the concentration of detected molecules) into an output one (the level of protein synthesised by the reporter plasmid). This processing device uses the continuous concentration of input metabolites as analog signals able to be processed by the metabolic or the genetic layer featured in the system. Multiple molecules can be continuously converted into the same effector at the same time allowing the construction of multi-input devices derived from analog adders.

Starting from this basis, we have implemented a machine learning architecture that allows integration of AI computations in cell-free systems [START_REF] Pandi | Metabolic perceptrons for neural computing in biological systems[END_REF]. One of the simplest architectures for signal processing is the perceptron [START_REF] Rosenblatt | The Perceptron, a Perceiving and Recognizing Automaton Project Para[END_REF]. A single-layer neural network transforms an input vector into an output value through the application of a weighted sum of the input followed by a thresholding function on the obtained result. Input data can be integrated as categorical or numerical (discrete or continuous). This signal-processing unit has been implemented in a cell-free system to enhance its multiplex sensing abilities [START_REF] Pandi | Metabolic perceptrons for neural computing in biological systems[END_REF]. This system is indeed capable of classifying complex samples in two categories (above or below a threshold of the reporter gene expression) based on the recognition of certain patterns of molecules present in those samples. A theoretical perceptron was first constructed and trained "in silico" to compute the optimal set of weights required to solve a defined problem of classification. The perceptron was then implemented in cell-free using a molecular model, giving concentrations of enzymes coding DNA that correspond to the weight of each input (see figure 1). The thresholding applied through the effector biosensor dose response curve.

[[ Figure 1 should be placed near here. ]] Different elements of the computing unit have a defined terminology: e.g. transducer, adder, actuator. A transducer, in cell-free terms, means one or more enzymes transforming one molecule into another. An adder consists of several reactions having different substrates but yielding the same compound. An actuator is simply a way for the cell-free system to report its activity, through the expression of a reporter gene, for example. Theoretically, many more parameters of a cell-free reaction could be implemented in such signal-processing manner.

The diversity of elements can be increased, as well as the number of layers in the perceptron, and the number of detectable molecules.

Here we describe a straightforward methodology that details how to use CAD tools to identify potential new biosensors for a chosen molecule, how to build and test these biosensor candidates and finally how to repurpose these cell-free biosensors into signal processing devices to implement neural computing in biological systems. This method to this date has been applied for the construction of a single layer 4-input perceptron but is easily scalable for more inputs and adaptable for more complex designs including multi-layer perceptron and other types of computational architectures. An example of potential implementation of a multi-layer perceptron using this system can be found in the supplementary figure 14 of the original metabolic perceptron paper [START_REF] Pandi | Metabolic perceptrons for neural computing in biological systems[END_REF].

Materials

Preparation of cell-free extract and buffer

Extract preparation equipment

In addition to common lab equipment, material for the preparation of Cell extract include large volume centrifugation equipment (for 1 L bottles) and a French Press.

Extract preparation media and buffers

Cells are grown in 2YTP media (31 g/L 2xYT, 40 mM potassium phosphate dibasic, 22 mM potassium phosphate monobasic) S30A and S30B buffers are also required for the preparation of the Cell Extract They are prepared with the following components:

 S30A = (

in silico materials

The R programming language (version 3.2.3) was used for fitting experimental data from actuator and transducers, as well as for defining and solving classification tasks. One can perform calculations described in section 3.6. with regular computational power, computational clusters are not required. Other programming languages can be used, but we recommend to use R in order to take advantage of our git repository containing all files from Pandi et al. work [START_REF] Pandi | Metabolic perceptrons for neural computing in biological systems[END_REF]. The repository is freely available here: https://github.com/brsynth/metabolic_perceptrons/tree/master/cell_free 3. Methods From this point using cross-references between databases or identifiers such as EC number, you can find enzyme candidates in UNIPROT with described catalysed reactions matching your expectations. You can also use the computational tool Selenzyme (http://selenzyme.synbiochem.co.uk/ ) with SENSIPATH provided references to find potential hits for other enzymes candidates (Note 3) 9. Retrieve DNA sequences for the identified parts. TF and Enzymes sequences can be codon optimised using any of the available tools. Promoter sequence is often defined as the 200 nucleotides before the start codon of a regulated gene. For Troubleshooting purposes, various size promoters can be synthesised to test for their response to the TF and their transcription initiation ability. (Note 4)

Constructing candidate biosensors plasmids

1. Synthesise the previously isolated sequence of your Transcription Factor, its regulated promoter and any enzyme required to convert the molecule you want to detect into the TF binding effector. This mix is then typically incubated in a thermocycler applying the following program: 30 s at 98 °C then 35 cycles with (10 s at 98 °C, 30 s at Tm +3 °C, 30 s/kb at 72 °C) then 2 minutes at 72 °C. Tm beings the lowest melting temperatures of the 2 primers and kb being the size of the amplicon in kilo bases. Use primer pairs FBP/RBP on the template plasmid 114598 to reamplify the linearized reporter backbone and the primers FBC/RBC on the template plasmid 114597 to reamplify the linearized backbone for enzyme or TF. Use the newly designed primers from step 2 on the synthesised parts from step 1 to reamplify the inserts. 4. Run an electrophoresis on the PCR product on a 1% Agarose gel stained with SYBR Safe using an appropriate DNA ladder to be able to discriminate your amplicon by its size. After an approximate time of 30 minutes at 100 V identify and cut the band of the expected size by imaging on a blue light transilluminator. 11. Finally spread 100µL of the final mix on LB Agar + Ampicillin (100 µL/mL) plates and incubate overnight at 37 °C.

12. Small scale culture for screening 1. Select 4 colonies per assembly and inoculate them in 3mL LB medium + Ampicillin (100 µL/mL) overnight.

2. Purify the plasmids from 2mL the bacterial cultures using a plasmid miniprep kit and use the primers seqF and seqR to check the integrity of each cassette by sanger sequencing.

3. Save one correct clone per construction by freezing the rest of the liquid culture at -80 ºC after addition of Glycerol to final concentration of 25%.

13. Large scale culture for plasmid production. You will require big quantities (>100 µg) of each plasmid to run the cell-free reactions for the characterisation and optimisation of each biosensor candidate (Note 7). To do so you have to realize large scale culture and plasmid extraction for each construction.

1. Inoculate 300 mL of LB + Ampicillin (100 µL/mL) from the -80 °C glycerol stock and grow the cells overnight.

2. Pellet the cells by spinning them at 6,000 g for 15 min at 4 °C.

3. Use the Maxiprep kit to recover plasmid DNA from the pellet. 4. After the last step of your purification, resuspend the precipitated DNA in 200 µL pure water in order to have a final solution at high concentration.

5. Measure the final concentration using a nanodrop and adjust it at 1 µM by either diluting it with water either concentrating it using a speedvac machine.

Preparing in house Cell-Free Extract and Buffer

The method briefly described here is adapted from a widely used protocol [START_REF] Sun | Protocols for Implementing an Escherichia coli Based TX-TL Cell-Free Expression System for Synthetic Biology[END_REF]of 3-PGA powered cell-free mix with minor modifications mostly concerning the lysis method and the starting strain.

2. Buffer preparation 1. Starting from individual solutions of chemicals dissolved in pure water (described in detail in the original paper [START_REF] Sun | Protocols for Implementing an Escherichia coli Based TX-TL Cell-Free Expression System for Synthetic Biology[END_REF]) prepare an Amino Acid (at 4X concentration) and an Energy solution mix (at 14X concentration) that will be used for buffer preparation. The Amino Acid mix has to be prepared by mixing all the 20 canonical Amino Acids at a final concentration of 6 mM except for leucine at 5 mM. Prepare Energy solution with HEPES pH 8 700 mM, ATP 21 mM, GTP 21 mM, CTP 12.6 mM, UTP 12.6 mM, tRNA 2.8 mg/ml, CoA 3.64 mM, NAD 4.62 mM, cAMP 10.5 mM, folinic Acid 0.95 mM, Spermidine 14 mM, 3-PGA 420 mM. Store each Mix at -80 ºC.

2. Using previously prepared extract, Amino acid Mix, and the Energy Solution, you will have to evaluate in this order the best concentration of: (i) Mgglutamate, (ii) K-Glutamate, (iii) DTT, and (iv) PEG8000 to add to your buffer to optimize the protein production of your extract. You will run 4 consecutives calibration experiments of 8h cell-free reactions. For each calibration, prepare a master mix for 12 reactions by mixing 88µl extract, 66µl Amino Acid mix, 18.86µl energy solution and 13.2 µL of the Addgene plasmid #40019 concentrated at 200 nM. Add the remaining three components that you are not calibrating for from the following four (Mg-glutamate, K-glutamate, PEG8000

and DTT) at either a starting concentration (if the optimum has not been determined yet), or the optimal concentration determined in a previous step. The concentration that leads to the highest GFP signal at 8h is identified as the optimal one to be used for future calibrations and run.

The starting concentration for each component are the following ones :(6mM for Mg-glu, 80mM for K-glu, 0 for DTT and 2% for PEG8000).

The tested concentration for each component are the following: Mg-glu : {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10 mM} K-glu : {0, 20, 40, 60, 80, 100, 120, 140, 160, 180 mM} DTT : {0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4 mM} PEG8000 {0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4%} 3. Prepare the final Buffer mix that will be used for running the cell-free reactions. The following quantities are for 1 tube of 650 µL (you need to scale up to have 2 tubes of buffer for 1 mL of final extract produced): per tube of buffer, add 128.97 µL of Amino Acid mix, 110.54 µL of Energy Solution, PEG8000, K-glu, Mg-glu and DTT according to the best determined concentration and pure water to adjust the volume to 650 µL.

Running Cell-Free reactions

All the cell-free experiments in the following parts should be run according to the same methodology. This method is for the preparation of N number of 20-µL cell free-reactions. All the following reactions have to be prepared on ice and have to be run in technical triplicate.

1. Thaw your cell free reagents (extract and buffer) on ice.

2. Pre-dilute your DNA and inducer stock: you will prepare 22-µL reactions each with 7.33 µL of extract, 9.17 µL of buffer and 5.5 µL of other components (DNA plasmids, inducer(s) and water to adjust). The easiest way to proceed if you have less than 5 different components to add per reaction is to make 20X solutions of the plasmids and inducers you want to use to add 1.1 µL of those to the reaction.

3. Prepare your cell-free Master Mix for (N+15%) reactions to compensate for the pipetting loss. Add (7.33*(N+15%)) µL of cell extract and (9.17*(N+15%)) µL of Buffer to a single 1.5-mL tube. 8. Pipette 20 µL from each tube to a well of a black 384-well plate pre-chilled. 9. Then cover the plate with a transparent sealing film.

10. Using a plate reader, monitor the green fluorescence (ex: 458 nm, em: 528 nm) at various gains during a 12-h kinetic run.

Cell-free biosensors characterizations and optimizations

In order to obtain the best possible biosensor response for our target we need to first develop an efficient biosensor for the TF binding effector, before optimizing the enzymatic conversion of our target molecule into our effector.

The first step is aimed to screen for any potential response of designed candidate biosensors and to answer 2 questions: first is there any interaction between the TF expression and the level of expression of the reporter and then is this interaction modulated in any way by the effector that we are trying to detect. To possibilities depending on the level of constraints you have with the input molecules. If you want to implement a computing device with a defined behavior but without any constraint on the nature of the inputs, you can reprogram the already developed benzoate based perceptron [START_REF] Pandi | Metabolic perceptrons for neural computing in biological systems[END_REF] to implement your desired computational function.

2. If you are interested in multiplexed sensing for defined molecules, you need to identify a common effector they can be converted into. You can use SENSIPATH with a reaction length of 2 to increase the chance to find a product accessible from the 2 or more substrates. Some central molecules (lactate, acetate, hydrogen peroxide, ammonium…) may be useful for that purpose but a specific attention has to be put on the development of biosensors with high dynamic range for these molecules with potential high noise issues related to their central positions on the metabolic networks.

  14 mM Mg glutamate, 60 mM K-glutamate, 50 mM Tris, pH 7.7) S30A is titrated using acetic acid, 2mM DTT are added just before use  S30B = (14 mM Mg glutamate, 60 mM K-glutamate, ~5 mM Tris, pH 8.2) S30B is titrated using 2 M Tris, 1 mM DTT are added just before use 3. Chemicals for the reaction buffer 1. Solutions of the 20 canonical amino acids, 19 of them concentrated at 168 mM, except leucine that is concentrated at 140 mM 2. Chemicals for the energy solution, including individual solutions of: HEPES pH 8, ATP, GTP CTP, UTP, tRNA, CoA, NAD, cAMP, folinic Acid, Spermidine and 3-PGA, all of them concentrated close to their limit of solubility.

3 .

 3 Additional components for the reaction buffer including K-glutamate solution at 3000 mM, Mg-glutamate solution at 100 mM, PEG8000 solution at 40% w/w and DTT solution at 1M. 4. Reference plasmid pBEST-OR2-OR1-Pr-UTR1-deGFP-T500 (Addgene #40019) for reaction buffer calibration experiments 2.2. Cloning of DNA parts and production of plasmids Additional materials for the construction and production of plasmid encoding the elements for Biosensors include : 1. High fidelity PCR Material (Q5 High-Fidelity 2X Master Mix or equivalent) 2. Template plasmid for backbone amplifications (plasmids #114597 and #114598 from Addgene) and respective primers : FWD primer Backbone for CDS : FBC: cccGGTCTCtGCTTactttatctgagaatagtc REV primer Backbone for CDS: RBC: cccGGTCTCtCATCatatctcttcttaaagttaaac FWD primer Backbone for Promoter : FBP:cccGGTCTCtATGCgtaaaggcgaagagctgttc REV primer Backbone for Promoter RBP:cccGGTCTCtTAAGaatagtaatacaggatccgaatcgtttcag 3. Thermocycler 4. TAE buffer, 50X 5. 1% agarose gel with SYBR™ Safe DNA Gel Stain or equivalent 6. Monarch® DNA Gel Extraction Kit or equivalent. 7. Thermo Scientific NanoDrop 1000 for DNA concentration determination or equivalent.

3. 1 .

 1 Identifying sensing routes for a target molecule [[ Figure 2 should be placed near here. ]] 1. Choose a target molecule to sense. Check that the molecule you want to sense isn't one of the components of the Cell-Free Buffer you are using. (Note 1) 2. Identify the InChi identifier or your target. To do so you can either search it by name on the database Pubchem ( https://pubchem.ncbi.nlm.nih.gov/ ) and search for the InChi in the section Identifier. The other possibility is to go on Pubchem sketcher (https://pubchem.ncbi.nlm.nih.gov//edit3/index.html) and to draw the molecule to retrieve it's InChI by replacing SMILEs by StdInChI. 3. Go to http://sensipath.micalis.fr/ to access the online platform. Use the section "Query with a Standard InChI" and paste the InChI of your molecule. 4. Run the tool first for one step pathways. If the results are not satisfying, you can run a second attempt for two step pathways. This can be the case if no pathway to a detectable molecule ( in green on the graph view) is identified or if the identified detectable molecule is non suitable for cell-free biosensors ( components of the Cell-Free Buffer) (Note 2) 5. Identify the promising effector molecules from the "Pathways view" section. If the molecule you want to sense appears in green, it means it is directly detectable without need for enzymatic conversion. For developing a biosensor for it you can skip the enzyme related steps. 6. Download the database of detectable molecules from: https://github.com/brsynth/detectable_metabolites and isolate the lines corresponding to the identified effector. First identify Transcription factor that interacts with your component. You can use column E to retrieve the names of it or, if this information is not available, use the column B that contains literature reference of papers describing potential sensing mechanisms for the molecule.7. Once you find a TF of interest you have to find the promoters that are potentially regulated by it. To do so you can search for it in databases such as RegulonDB (http://regulondb.ccg.unam.mx/) for E. coli, Subtiwiki (http://subtiwiki.unigoettingen.de/) for B. subtilis. You can also try a naive bibliographic search for this regulator to identify features linked to it like regulated promoters but also mechanism of action and possible existing design of biosensors using it. 8. Identify the enzymes converting the molecule you want to sense into the effector for which you found a TF. By clicking on the edges from the SENSIPATH graph view you can retrieve the references associated with the enzymatic reaction of interest.

2 .

 2 Design and orders primers for golden gate cloning of the synthesised parts. Genes (TF and enzymes) requires overhangs in the format, FWD : ccGGTCTCtGATG…. REV : ccGGTCTCtAAGC…. promoters requires overhangs in the format, FWD : ccGGTCTCtCTTA…. REV : ccGGTCTCtGCAT…. 3. Run High fidelity PCRs to amplify with the correct Golden Gate overhangs the Vectors and the synthesised inserts. The reaction typically consists of pipetting 25 µL of 2X Q5® polymerase Master Mix, 2.5 µL of FW and RV primer at 10 µM, 1µL of template DNA (at concentration around 100ng/µL) and 19 µL of water.

5 .

 5 Recover the DNA from the Gel using a DNA Gel extraction kit, following the kit's instructions to purify your DNA fragment. Determine the titer of purified DNA using a nanodrop spectrophotometer. 6. Prepare a golden gate reaction to insert each fragment in its respective backbone (Note 5). To do so you need to calculate the molarity of your DNA. You can use the tool available at https://nebiocalculator.neb.com/#!/dsdnaamt using the nanodrop determined concentrations of each purified DNA fragment. For the golden gate reaction incubate 100 fmol of insert with 50 fmol of linearized backbone in a tube with 1 µL T4 DNA ligase, 1 µL BsaI enzyme, 2 µL T4 DNA ligase Buffer and water to adjust the volume to 20 µL. Incubate the mix at the 37 °C for 1 h and then 16 °C for 5 minutes. (Note 6) 7. Transform DH5α competent cells with the Golden Gate Reaction product. Incubate 5 µL of golden gate product with 50 µL chemically competent cells at 4° C for 30 min. 8. Heatshock at 42 °C for 45 s. 9. Incubate at 4 °C for 3 min. 10. Add 300 µL LB media and incubate at 37 °C for 1h.

Finally, add water

  to this mix to reach a final volume of 237.6 µl. Prepare each calibration reaction by mixing in PCR tubes 19.8 µl of the master mix with 2.2 µl of the tested component concentrated at 20X. Then, take 20 µl from each PCR tube and pipette them inside individual wells of a 384 wells microplate before incubating the plate for 8h while measuring GFP signal produced (ex: 458 nm, em: 528 nm).

4 .

 4 Mix by briefly vortexing. 5. Pipette 16.5 µL of that mix in N PCR tubes (in strips).6. Add the respective other inputs (plasmids, inducers, water…) from the 20X stocks, to each PCR tube up to 22 µL. 7. Close these tubes and mix them by briefly vortexing and bench centrifugation.

2 .

 2 answer those, you need to characterize the behavior of the reporter plasmid in a cellfree reaction in presence or absence of Transcription Factor and in presence or absence of Effector. If you have identified multiple candidates for the sensing of the same molecule this first experiment can also be used to decide what are the TFs and the promoters most promising for the nexts steps. 1. Run cell-free reactions varying the quantity of added plasmids and chemicals in the following possibilities: Reporter DNA concentration at 0 or 30 nM, TF DNA concentration at 0 or 30 nM and inducer concentration at 0, 100 µM or 1 mM. Use end point kinetics results at 8 hours to evaluate the potential mechanism of your TF (activation or repression) and a possible response of your system to the inducer. 3. Find the optimal plasmid concentration for the Reporter and TF that results in a maximal response for your biosensor. With 3 concentrations of inducer (0 µM, 100 µM and 1 mM), test a combination of concentrations gradients for the added plasmid DNA varying reporter and TF DNA concentration on a logarithmic scale (0, 0.1, 0.3, 1, 3, 10, 30 and 100 nM) 4. Using the previously determined best pair of plasmid concentration, evaluate the inducer dose response of the constructed sensor. Run a cell-free experiment with inducers concentration at : 0 nM, 0.5 nM, 1 nM, 2 nM, 5 nM, 10 nM, 20 nM, 50 nM, 100 nM, 200 nM, 500 nM, and 1,000 nM. Dose response curve is obtained by plotting (fluorescence at concentration x / fluorescence at concentration 0) for each datapoint. Use the final results obtained to assess the performance of the developed sensor and its conformity with your objective of sensitivity. (Note 10) 5. Once you have a satisfying sensor for your effector, you can start screening for enzymes candidates to convert your target into this effector. The first experiment to do consists of screening for any potential activity on each of the selected enzymes candidates. Prepare a cell-free experiment with the previously identified best concentration for reporter and TF DNA, varying Enzyme DNA concentration at 0 nM or 10 nM and the molecule of interest at 0 mM or 1 mM. 6. Select the best enzyme candidate and optimize the expression of the enzyme on a dose response curve (same concentrations as in step 3) by varying the DNA concentration coding for the enzyme on the scale of 0, 0.1, 0.3, 1, 3, 10, 30, 100 nM.

3. 6 .

 6 Design, Build and Test a perceptron 1. Define the architecture of the perceptron based on input molecules. There are 2

3 .

 3 Create and fit the cell-free model. It can be decomposed in 2 parts: actuators and transducers.4. An actuator is modelled with a modified Hill function; commonly used in biochemistry for simulating the binding of ligands to proteins according to the ligand concentration.

-

  total: concentration of input metabolite in µM -: concentration of input metabolite yielding half of the maximum induction of the system (also called IC50) -ℎ : Hill coefficient characterizing the cooperativity of the induction system -fc : stands for fold-change, corresponds to the dynamic range (in Arbitrary Units) of the system -basal : basal GFP fluorescence signal without input, i.e. the background noise of the system -lin : accounts for the linearity of the system when dealing with concentrations saturating the Hill transfer function

Acknowledgments

PS and JLF are supported by the ANR SynBioDiag grant number ANR-18-CE33-0015. LF is supported by the French National Research Institute for Agriculture, Food, and Environment (INRAE), through the "Métaprogramme BIOLPREDICT". MK acknowledges funding support from Ile-de-France region's DIM-RFSI, INRAE's MICA department and the University of Paris-Saclay.

We recommend at least some inspiration from our fitting process described in the notes section (Note 11) as it puts an emphasis on tackling the loss of signal of the actuator when the whole system is modelled then implemented.

Otherwise, one can adapt another fitting process for a particular project's needs.

6. Measure the model's performance with different metrics: Root Mean Square Deviation (RMSD); R²; Weighted R² ; Error percentage...

Ensuring high metrics (e.g. above 0.9 for scores between 0-1) on many experimental data points guarantees a robust model for predicting the weights schemes for each classification task.

7. Define your classification problem. Start by defining a set of tasks (here, classifications), which outputs either 0 (OFF) or 1 (ON); as well as a set of weights to be tested (equivalent to a range of possible enzyme concentrations in the cell-free experiment, for our work it was between 0.1-10nM).

8. Continue by sampling uniformly input values for your problem, that is to be resolved by the perceptron. For example, in the case of a binary classification of solutions composed of hippurate and cocaine, sample points in a given range, either for a "low" concentration or a "high" concentration for each of the compounds. Here, let us assume we can sample between 0 and 2 µM for low concentrations and between 80 and 100 µM for high concentrations.

9. For each of the two clusters, choose to sample either in the low or high range for each compound.

10. Then, two clusters have been produced and a binary classifier can be easily defined on these points: the perceptron set of weights and its corresponding fluorescence threshold. Please find a visual example of the definition of a classification task on figure 3. Further is detailed how to find the set of weights and the fluorescence "decision threshold".

[[ Figure 3 should be placed near here. ]]

11. Predict the best set of weights for solving this problem (aka "train" the perceptron on a classification task). Using the previously fitted (part 2) and benchmarked (part 3) model, simulate all possible input combinations with all possible sets of weights.

12. Screen for performant set of weights, i.e. those enabling a sharp threshold between the output states of the system ("ON" or "OFF"). Several thresholds can be tested for considering an output value as ON or OFF, for all possible simulations.

13. Manually select the best set of weights and corresponding threshold, i.e. those that show the highest and clearest difference between ON and OFF behaviors, and in most scenarios. Also, prefer those showing low enzyme concentrations (to avoid resource competition). One can also test several possibilities (several sets of weights and/or thresholds) in the following cell-free experiments.

14. Implement the designs predicted as best ones, to test it in a cell-free reaction.

Start by drawing a test set of chemical (input) combinations from the compositions used for the training of the in silico perceptron. Be cautious to evaluate enough points from the space of possibilities to capture the behavior of your system. If your perceptron is designed to solve the particular problem of binary inputs sample classifications, you can eventually build a complete test set with all the possible combinations of inputs (see figure 1D).

15. The designed test set can be prepared at 20X concentration in PCR tube strips to be evaluated on various perceptron implementations. We advise you to prepare a single master mix with the cell-extract, the buffer and the DNA coding for the reporter system and the various transducers at the desired concentration for each implementation of the perceptron that you designed. [START_REF] Pandi | Metabolic perceptrons for neural computing in biological systems[END_REF]. Use this master mix to evaluate the response of the system on the chemical test set previously designed. If the perceptron does not have the expected behaviour from in silico analysis, try another set of weights.

Notes

1. Most of the molecules present in the cell-free buffer are in high concentration masking any potential response for a biosensor designed to detect them. This list includes all the 20 amino acids, HEPES, ATP, GTP, UTP, CTP, tRNA, CoA, NAD, cAMP, folinate, spermidine, 3-PGA, Magnesium, Potassium and DTT.

2. If you cannot find satisfying results with the software SENSIPATH you can try to find other SEMP using the retrosynthesis workflow Retropath [START_REF] Delépine | RetroPath2.0: A retrosynthesis workflow for metabolic engineers[END_REF]. It has the advantage to allow prediction for pathways with more than 2 steps or pathways using promiscuous activity of enzymes to find potential new reactions. To run retropath for biosensor design, use the list of detectable molecules [10] as Sink, the molecule you want to detect as source and the reaction rules in the forward direction. Results coming from this workflow have to be taken with more care as it relies on less reliable predictions.

3. The Selenzyme tool is predicting potential enzymes catalysing a defined reaction based on similarity of sequences, reactions or other features existing between enzymes from a well-annotated database and enzymes potentially catalysing the query reaction.

The predictions should be manually checked in published literature/ datasets to evaluate if the identified enzyme can likely catalyse the given reaction.

4. The methodology described here is a standard to be used in the case where the identified promoter has limited annotation or features described. You are strongly encouraged to search for existing biosensing projects described in the literature that uses the same Transcription Factor with defined size promoters or synthetic ones (built by inserting TF binding sequences in another promoter) as these promoters may show a better response than the natural ones.

Golden Gate Assembly was chosen over other methods like Gibson Assembly as it

allows reusing the same primers to reamplify backbone for every new insert cloned.

6. This thermocycler protocol for Golden Gate is a variation of the fast Golden gate assembly protocol (1h 37 °C, 10 min 55 °C) that is adapted to backbone plasmid containing internal bsaI cut sites. Removing the last 55 °C step and adding one at 16

°C for 5 min avoid cutting the final assembly containing the BsaI site. First, fit the actuator experimental data (effector dose response curve obtained in part 3.5 section 3) to the Hill function model.

To do so, fit 100 times the actuator experimental values to the actuator model. Let all parameters be able to vary. Use ordinary least-square error or the R "optim" function (that uses the Limited-memory Broyden Fletcher Goldfarb Shanno algorithm) as the objective function for the data fitting. Keep the seed of the fitting process in memory so it ensures reproducibility. Retrieve fitted parameters for each fitting (a population of 100 sets of fitted parameters will be produced). From this population, save the mean, standard deviation, standard error and confidence interval; for each parameter.

To account for the decrease of signal in experimental actuator data when the whole system is implemented in cell-free, one needs to fit the transducer model in a specific way.

To fit the transducers models, we need to actually fit the whole model to the whole system's cell-free experimental data.

In order to have a coherent perceptron, we will actually constrain the actuator's parameters, then fit all transducers and the actuator together.

To do so, start by constraining each of the actuator's previously fitted parameters with bounds corresponding to the 95% confidence interval; or +/-one standard deviation from the mean, in the case of following parameters: the fold-change, termed "fc", and the baseline, termed "basal".

We call "pseudo-experimental" the data used to fit the whole system, but that originates from specific components' experimental data. These individual data were aggregated as described in the whole-system model: combine transducers data by simply adding their output and feeding the sum to the actuator. Note that in the cellfree experimental system, the same process happens: e.g. benzoate concentrations are added as a result of all transducers yielding benzoate, and the actuator takes this aggregated benzoate concentration as an input. As a result, we obtain "pseudoexperimental" data for the whole cell-free system, with experimental data only available for each component separately.

Once your whole system is modelled and you have constrained your actuator's parameters, you can fit your transducer parameters. Initialize the whole model by drawing actuator parameters values according to a Gaussian distribution centered on the mean of the parameters estimations, with a standard deviation equal to the standard error of this parameter estimation. Let all transducer parameters vary to fit the whole system to the previously described 'pseudo-experimental' data.