## 2 Materials and Methods  
### A General Framework for Individual Subsystem Models  
General ISMs give rise to a system of implicit Differential Algebraic Equations (DAEs) for each subsystem [47] of the form:  
\[\Phi_{i}\big{(}t,x_{i}(t),\dot{x}_{i}(t),y_{i}(t),v_{i}(t),u_{i}(t),p_{i}\big{)} =0, \tag{1}\]  
where \(t\) is time, \(x_{i}\) are the differential variables, \(y_{i}\) are the subsystem output algebraic variables, \(v_{i}(t)\) is the subsystem input vector (that is, inputs whose values are equal to outputs of other subsystems in the ISM), \(u_{i}\) is the ISM external inputs vector (that is, the inputs, if any, coming from the outside environment), and \(p_{i}\) is the subsystem internal constant parameters vector. In DAEs, some of the equations could be purely algebraic, leading to more difficulties when combined with differential equations [48]. The interconnection of subsystems is achieved by connecting subsystem outputs to the inputs of other subsystems. Precisely, the \(i\)-th subsystem is connected to the \(j\)-th subsystem if there is a directed connection from one of the scalar components \(\widetilde{y}_{ir}\) of the output vector \(y_{i}=(\widetilde{y}_{ir})\) of the \(i\)-th subsystem to a scalar component \(\widetilde{u}_{js}\) of the input vector \(u_{j}=\big{(}\widetilde{u}_{js}\big{)}\) of the \(j\)-th subsystem. Without loss of generality, we can neglect self-loops (i.e., we neglect the case \(j=i\)) using the identity function as an intermediate auxiliary subsystem. Appending vectors \(x_{i}X=[x_{i}]^{\widetilde{\text{\tiny{i}}}\widetilde{\text{\tiny{i}}} \widetilde{\text{\tiny{i}}}}\) and the same can be done \(Y=[y_{i}]\), \(V=[v_{i}]\), \(U=[u_{i}]\), \(P=[p_{i}]^{\widetilde{\text{\tiny{i}}}\widetilde{\text{\tiny{i}}}}\)\(\Phi=[\Phi_{i}]^{\widetilde{\text{\tiny{i}}}\widetilde{\text{\tiny{i}}} \widetilde{\text{\tiny{i}}}}\). Therefore, combining all the subsystem equations, the result is a fully implicit system of DAEs that is generally solved with different available methods, including symbolic manipulation and numerical solvers:  
\[\Phi(t,X,\dot{X},Y,V,U,P)=0. \tag{2}\]  
On the other hand, without loss of generality, we can describe our method using the Ordinary Differential Equations ODE representation for better readability and using a suitable graphical representation. DAEs can be easily introduced without any change in the method. Moreover, if the Jacobian of \(\Phi\) with respect to \(\dot{X}\) is non-singular, all the following forms are equivalent to (2) [47]. In particular, for each subsystem \(N_{i}\) we have:  
\[f_{i}\big{[}t,x_{i},\dot{x}_{i},v_{i},u_{i},p_{i}\big{]}=0, \tag{3a}\] \[y_{i}(t)=g_{i}\big{[}t,x_{i},\dot{x}_{i},v_{i}\big{]}. \tag{3b}\]  
All the algebraic variables disappear and the \(x_{i}\) are all differential ones. Sometimes ODEs can be written in the explicit form:  
\[\dot{x}_{i}=f_{i}(t,x_{i}(t),v_{i}(t),u_{i}(t),p_{i}). \tag{4}\]  
Such an interconnected ISM formed by all the subsystems \(N_{i}\) can be represented as shown in Figure 1. As said above, appending vectors \(v_{i}=[v_{ir}]\) allows to form a single column vector \(V\):  
\[V(t)=[v_{1}(t),v_{2}(t),\dots,v_{N}(t)]^{T}=[V_{h}], \tag{5}\]  
where \(V_{h}\) is the \(h-\)th scalar component of \(V\); the same can be conducted with outputs \(y_{i}=[y_{is}]\):  
\[V=MY \tag{6}\]where each element \(M_{hk}\) of the connection matrix, \(M\) is \(M_{hk}=0\) or \(M_{hk}=1\). Last condition, \(M_{hk}=1\) represents a connection from the scalar output \(Y_{k}\) to the internal scalar input \(V_{h}\). If no self-loops are allowed, all the diagonal elements are null: this hypothesis can be accepted as it simplifies the calculations below without losing generality. Finally, we can assume that each scalar component of a subsystem output vector is connected to the one and only one scalar component of another subsystem internal input vector (multiple connections can be modelled with multiple identical scalar outputs). To this end, it helps in explaining the new approach in the next section.