Mathematical models contribute to our understanding of Biology in several ways ranging from the 2 quantification of biological processes to reconciling conflicting experiments [1] . In many cases, this requires 3 formulating a mathematical model and extracting quantitative estimates of its parameters from the 4 experimental data. Parameters are typically unknown constants that change the behavior of the model. 5 While it is usually recognized that parameter estimation requires the availability of sufficient informative 6 data, sometimes it is not possible to estimate all parameters due to the structure of the model (whatever the 7 quantity or quality of the data). This inability is referred to as 'structural identifiability', a concept 8 introduced decades ago by Bellman andÅström [2, 3] , as opposed to the 'practical identifiability' that 9 depends on limitations set by the data. Practical identifiability has important consequences that can lead to 10 questionable interpretations of the data leading to some recent controversy around this point [4, 5] . 11 Structural identifiability is a necessary condition for model fitting and should be used before any attempt to 12 extract information about the parameters, and as a test of the applicability of the model itself. Structural 13 1/10 identifiability can be qualified as global or local. Global structural identifiability tests the ability to estimate 14 unique sets of parameters, while local (or simply, structural identifiability) means that parameters can be 15 estimated only in a limited subset of the space of parameters, i.e., only combinations of parameters are 16 identifiable [6] [7] [8] [9] [10] . In practical terms, these definitions can be translated into the language of sensitivity 17 analysis as identifiability requires that (i) the columns of the sensitivity (or, equivalently, the elasticity) 18 matrix are linearly independent, and (ii) each of its columns has at least one large entry [11, 12] . 19 Traditionally, work primarily focused on linear systems [2, 3, 13] based on ordinary differential equations 20 (ODE). For non-linear models, those methods cannot be applied, so many methods have been proposed in the 21 literature to address structural identifiability. Early attempts were based on power series expansions of the 22 original non-linear system [14] , the similarity transformation method [15] [16] [17] or the so-called direct-test 23 method proposed by Denis-Vidal and Joly-Blanchard [18, 19] . These methods exploit the definition of 24 identifiability either analytically [18] or numerically [20] [21] [22] [23] [24] [25] , but they are not generically suitable for 25 high-dimensional problems. Xia and Moog [6, 26] proposed an alternative to these classical methods based on 26 the implicit function theorem, but this method also becomes involved to apply for complex models [27] . 27 Another approach that is becoming mainstream is based on the framework of differential algebra [28] [29] [30] [31] . 28 These methods are also difficult to apply, requiring advanced mathematical skills and, in some cases, replace 29 highly non-linear terms by polynomial approximations that simplify the analysis. On the positive side, they 30 are based on rigorous mathematical theories, are suitable for non-linear models and, more importantly, they 31 can be coded using existing symbolic computational libraries. In this regard, it is worth mentioning 32 DAISY [32] , GenSSI [33] , COMBOS [34] or, more recently, SIAN [35] . 33 In almost all cases, the major disadvantage of these methods is their difficulty to apply them to even a few 34 differential equations, hence requiring advanced mathematical skills and/or dedicated numerical or symbolic 35 software (that is frequently unable to handle the complexity of the problem). This explains why, despite the 36 huge volume of publications in the field of theoretical biology, only a few address parameter identifiability 37 explicitly. In this paper, we introduce a simple method to assess local structural identifiability of ODE 38 models that reduces the complexity of existing methods and can bring identifiability testing to a broader 39 audience. Our method is based on simple scaling transformations, and the solution of simple sparse systems 40 of equations. Identifiability for stochastic models [36] is out of the scope of our work. Consider a simple death model in which the death rate is the product of two parameters λ 1 and λ 2 , namely
with the solution
It is evident that from an experiment only the product λ 1 λ 2 can be inferred, and not any of the two 46 independently. Following the 'actionable' definition in Ref. [11] , local structural identifiability is directly 47 linked to the linear independence of the columns of the elasticity (or, sensitivity) matrix. In this case, the 48 elasticity matrix would be simply a 1 × 2 matrix,
We now propose to multiply λ 1 with a generic scale factor u, and to divide λ 2 by the same factor, such that 51 the solution remains invariant. Deriving the scaled solution of eq. (2) with respect to that scale factor u, and 52 by the chain rule, 53 dx du = 0 (as u is arbitrary) (4) and, also,
where the last equality follows from Eq. (4)
Rearranging Eq. (5) and dividing by x,
so both columns of the elasticity matrix are linearly dependent and, accordingly, λ 1 and λ 2 are unidentifiable. 57 In this case we had complete knowledge of the solution, and consequently, it was straightforward to find the 58 right way to introduce the scaling u. Fortunately, this simple scaling calculation can also be performed 59 directly on eq. (1). Introducing two unknown scaling factors, u 1 and u 2 , into that equation,
Requiring that this remains identical (or, more formally, invariant) to Eq. (1), i.e., λ 1 λ 2 x = u 1 λ 1 u 2 λ 2 x, 61 hence u 1 u 2 = 1. The fact that u 1 and u 2 cannot be solved individually, also means that the real values of λ 1 62 and λ 2 cannot be determined, namely both parameters are unidentifiable.
Next consider a death model with immigration:
In this case, to leave the system invariant we need to find u 1 and u 2 such that
for all values of x at any time. Rearranging the latter equation,
where the left-hand side of the last equation is a constant and the right-hand side depends on time. Hence the 67 only possible solution to the latter equation is u 1 = u 2 = 1 implying that both λ 1 and λ 2 are locally 
where the functions f i depend on the specific details of the problem at hand and x i,0 are the initial 75 conditions. We need to distinguish between those variables that can be observed (measured) in the 76 experiment, x 1 . . . x r , and those which cannot (they are often referred to as latent variables), x r+1 . . . x n .
According to Sec. 1 in the Supporting Information, each function f i is split into M functional independent 78 summands, f ik ,
having the property that f ik is functionally independent of f il for every k = l. 1. Scale all parameters and all unobserved variables by unknown scaling factors, u:
and substitute them into Eqs. (10) below.
2. Equate each functionally independent function, f ik , to its scaled version. Namely,
where u xi = 1 for 1 ≤ i ≤ r and the prefactor in the right-hand side of the equation comes from the scaling of dxi dt → u xi dxi dt . From Eq. (8b) it follows that u xi = u xi,0 . 3. From the Eqs. (9) , find combinations of the scaling factors u that leave the system invariant.
Hereafter, we will denote these as the identifiability equations of the model. We summarize our method in Box 1.
In summary, our method reduces the complexity of finding identifiable parameters to finding which scaling 86 factors do not satisfy the trivial solution u i = 1. In the literature, when a scaling factor is related to one of 87 the latent variables x r+1 . . . x n , if u x k = 1, then x k is said to be observable [10] . Thus, our method addresses 88 at the same time identifiability and observability. Additionally, irreducible equations involving two or more 89 parameters provide the so-called identifiable groups of variables that cannot be fitted independently. In the 90 case of the pure death model above, the identifiability equation u λ1 u λ2 = 1 is a signature of the unidentifiable 91 group λ 1 λ 2 . This is interesting as groups involving latent variables (for instance, u xj u λ k ) would inform future 92 experiments aimed to observe that variable and decouple that group.
It is also worth mentioning that our identifiability test (illustrated by example in the Supporting Information) 94 provides a simple way to find a type of symmetry that is related to scale invariance. More sophisticated 95 methods have been introduced in the literature to address other symmetries [37] [38] [39] using the theory of Lie 96 group transformations, however, that approach involves complex calculations assisted by symbolic 97 computations.
The main result 100 Consider a model described by a set of n ordinary differential equations (ODE)
where f ik is functionally independent of f il for every k = l (namely, they satisfy the generalized Wronskian (1)-(3), we seek for scaling of the parameters that leave the system invariant. As we prove 105 below, this invariance (or lack of) is related to the identifiability of the parameters. Hence, if we define the 106 following scaling transformation:
(where the variables x 1 . . . x r are unmodified as we can measure them in the experiment) we can write the following set of re-scaled equations:
where M is the number of functional independent summands in the equation. It is convenient to rewrite 108
Eq. (12c) as
to perform the scale invariance analysis below in a simpler way.
If the solution is invariant under this transformation, then the right-hand sides of Eq. (10) and, consequently 111 Eqs (12) should be equal. Besides, by the functional linear independence of the functions f ik we can split 112 each summand. Thus,
and
These new set of equations are much easier to solve than the ones that we would obtain from Eqs (12a)-(12c) 115 (which would be equivalent to the so-called direct-test method [18] ). 116 We can express the solution of these equations as 117 u λ k = F (u xm 1 , u xm 2 , . . . , u λj 1 , u λj 2 , . . .),
We denote these the identifiability equations of the model. For each parameter k, the identifiability 118 equation will depend only on a few other scaling factors m 1 , m 2 , . . ..
If take the partial derivative of the transformed solution 120
x i (x 1 , . . . , x r , u xr+1 x r+1 . . . u xn x n ; u λ1 λ 1 , . . . , u λm λ m ) with respect to u λ k , we find (by the chain rule) 
At this point we have two possibilities. Either all the coefficients β kj = 0, meaning that Eq. (16) is not 127 satisfied (we will simply have u λ k = 1), or the columns of the elasticity matrix related to the parameter λ k 128 has linearly dependent columns, and then the parameters λ k , λ j1 , λ j2 , . . . are unidentifiable.
On the other hand, if the solution of the identifiability equation is u λ k = 1 for some parameter, then the 130 column of the elasticity matrix coming from that parameter is linearly independent of the others, and the 131 parameter is local structurally identifiable. The adjective "local" is required because the method stems on 132 the continuity of the derivative of the solution of Eq. (17) . Thus, it is unable to capture any discrete 133 (countable) transformations like, for instance, those related to exchanging two parameters of the model.
Finally, by Eq. (12d) we find that u xi = u xi,0 , so an initial condition is identifiable whenever its 135 corresponding variable is observable.
Comparison with other methods 137 We have applied the method outlined in Box 1 to 13 different models defined and analyzed in detail in the 138 Supporting Information. The choice is based on two criteria: on the one hand, models 1-5 are included for 139 pedagogical purposes. They are simple enough to illustrate the method and most of the existing methods also 140 provide definite answers. Models 6-13 were chosen because they have previously been analyzed using the 141 methods summarized in the Introduction and in Table 1 . This allows us to put our method in direct 142 competition with those methods and to highlight their merits and limitations.
The results of this comparison are summarized in Table 2 , which is inspired by a similar table in Ref. [7] .
144 Table 1 . List of current methods testing structural identifiability. We introduce here the acronyms referred to in Table 2 . Information, and that, after some practice (and using some interesting motifs as having sums of different 153 parameters, or the coefficients related to diagonal terms in the system of equations) the calculations can be 154 Table 2 . Summary of models compared in the literature: The number in brackets in the Model Name column correspond to the number of observed variables. Model Numbers correspond to those in Table 1 in the Supporting Information. The acronyms for the methods are summarized in Table 1. This table has been  inspired by Table 1 in Ref. [7] . [27] or take hours using symbolic computation packages.
Together, this broad applicable and simplicity are the main features of our method and this may attract the 157 interest of mathematical modelers and spread the culture of checking structural identifiability as a mandatory 158 step when fitting experimental data. 159 We would like to highlight a connection with the so-called Buckingham-Π theorem of dimensional 160 analysis [48] . In some sense, the scale invariance property is related to the principle of dimensional 161 homogeneity, i.e., the constraints on the functional form of the independent variables with the parameters.
Our identifiability equations are therefore similar to finding the so-called Π-groups in the theorem.
A limitation of the method is that it is restricted to testing local identifiability. This is implicit in the 164 differentiability of the elasticity matrix which, by definition, is a local operation. Discrete symmetries are not 165 captured, and more sophisticated methods (based on Lie group transformations [39] ) are required. However, 166 simple manipulation of the equations to remove the latent variables can improve the explanatory power of the 167 method and might capture those discrete symmetries (see Sec. 3.8 of the Supporting Information). We leave 168 that extension for future developments.
Finally, in this work we have chosen to solve the scaling factor equations directly as it is easy to perform with 170 pen and paper. However, if we were to redefine the scaling factors as u i = e wi , the new factors w i would obey 171 a linear system of homogeneous equations. It is therefore expected that the problem of identifiability is related 172 to the rank of the matrix defining the linear system of equations. In that regard, the theorems presented in 173 the Supporting Information could be supplemented with generic results on homogeneous systems of equations. 174 Thus, our results provide a solid ground for the method and indicate a venue for further development in other 175 systems like delay-differential or partial differential equations. Finally, while we emphasize the simplicity of 176 the method, it is obviously amenable to be implemented using symbolic computation packages.
In the Supporting Information we collect the theorems sustaining the method and a catalogue of models with 179 a detailed computation of the identifiability equations that were used to build Table 2 . 