Cash Flow Optimization for Renewable Energy Construction Projects with a New Approach to Critical Chain Scheduling

: This study concerns the use of the critical chain method to schedule the construction of renewable energy facilities. The critical chain method is recognized as a useful project management tool, transforming a stochastic problem of uncertainty in activity durations into a deterministic one. However, this method has some shortcomings. There are no clear principles of grouping non-critical activities into feeding chains. Another ambiguity is sizing the feeding buffers with regard to the topology of the network model and the resulting dependencies between activities, located in different chains. As a result, it is often necessary to arbitrarily adjust the calculated sizes of feeding buffers before inserting them into the schedule. The authors present the new approach to sizing the time buffers in the schedule, enabling a quick assessment of the quality of a given solution variant and ﬁnding a solution that best meets the established criteria, conditions, and constraints. The essence of the presented approach is the two-step sizing of time buffers with the use of deterministic optimization and stochastic optimization techniques. Taking into account construction management needs, the optimization criteria are based on the construction project cash ﬂow analysis. The effectiveness of the presented approach is illustrated by an example of developing a wind power plant construction schedule. According to the results, the presented approach ensures the protection of the scheduled completion date of the construction and the stability of the schedule.


Introduction
Wind is a renewable energy source that is completely free and inexhaustible. Moreover, wind energy is widely available and therefore does not need to be imported. The process of generating such energy does not generate side effects in the form of harmful carbon dioxide emissions and other pollutants to the atmosphere. Wind farms have a high energy potential and are able to generate up to 70 times more energy than is needed for their installation, operation, and dismantling. The process of building such farms can, however, be optimized, which would allow for even greater profitability of such investments and counteracting disruptions during their implementation. This topic is particularly important from the point of view of the increasing use of energy sources and progressing climate change [1][2][3].
However, the investment optimization process is complex [4]. Construction projects are characterized by the uncertainty of implementation conditions, despite the repetition of individual components of the investment process [5][6][7]. One of the consequences of this uncertainty is the incorrect estimation of various disturbances, causing delays in the performance of individual works and the entire project, and the related financial losses of the contractor. These losses are mainly caused by contractual penalties charged by the ordering party and subcontractors (if applicable).
Additionally, in construction planning practice, the contractor often does not have sufficient renewable or non-renewable resources to complete the works in the shortest possible time [8]. This may be due to limitations in the availability of workers with specific skills, limitations in the availability of certain units of construction equipment, or limitations in the availability of capital to cover current construction expenses. Limitations in the availability of resources may apply to the entire construction period or to its individual subperiods. Limitations may also apply to the level of resource consumption. For example, the limited space of the construction site may decrease the allowable daily number of workers employed on the construction site due to the limited space allocated to social and administrative temporary buildings on site. Therefore, it often becomes necessary to link the periods of execution of the individual works with the availability of resources used for their execution [9].
The critical chain method is recognized as a useful project management tool, capable of dealing with such problems. However, this method has some shortcomings, which are addressed in this manuscript.

The Conventional Approach to the Critical Chain Method
Numerous literature sources point out that, when taking into account resources constraints and uncertainty of activities' durations, the critical chain method (CCM) is a proper tool for ensuring efficient resource management and stability of the project completion time [10]. In the CCM, the stability of the project completion date is achieved by the following: • The identification of activities that form the critical chain and feeding chains; • The insertion of two specific kinds of time buffers (safety margins) into the schedulethe project buffer, stabilizing the completion date of the entire project, and feeding buffers, protecting the critical chain against the propagation of disturbances through the individual feeding chains.
The critical chain may be defined as a resource-constrained critical path in the network model, i.e., as a set of activities, whose sum of durations determines the duration of the entire project, taking into account technological sequence of activities and additional, informal relations introduced in the network model to balance renewable resources with limited availability (resource relations). The starting dates of activities in the critical chain are scheduled as soon as possible. Feeding chains group non-critical activities. Their starting dates are scheduled as late as possible.
There are two stages in scheduling a project using the critical chain method. In the first stage, the initial schedule is prepared. For the preparation of the initial schedule, "pessimistic" estimates of activities' durations are used, considered to be safe. Using the critical path analysis, the latest times for the execution of individual activities are determined. Then, the renewable resources are leveled, taking into account the limitation of the availability of renewable resources. At the end of the first stage, the critical chain and feeding chains are identified. In the second stage, firstly, the scheduled durations of works are reduced from "pessimistic" to "aggressive" values.
In numerous literature sources, the median of the probability distribution of activity duration is suggested as an "aggressive" estimate. However, in practical applications, half of the "pessimistic" estimate is often used as an "aggressive" estimate, according to the original Goldratt's proposal. After reducing the activities' durations, a buffered project schedule is prepared. The project buffer and feeding buffers are calculated by aggregating the time reserves of activities forming the critical chain and feeding chains.
Most literature sources suggest two analytical methods of calculating the sizes of buffers: the cut-and-paste method and the root-square-error method. There are other, numerous approaches to the buffers' sizing. Some of them are based on the number of precedence relationships, resources tightness (the closer the resource usage is to the resource availability, the more likely a delay will occur) [11] risk class assessment, or expected and pessimistic activity durations. An interesting research study that compares buffer-sizing techniques in the CCM is presented in [12] with the conclusion that the best method for buffer sizing is stochastic optimization.
It should be noted that the CCM focuses on the stability of the completion date of the entire project. In order to ensure the stability of starting dates of individual activities, the method presented by [13] or by [14] can be used. The essence of both methods is the appropriate sequencing of activities in the schedule, taking into account the limited availability of renewable resources. However, these methods are differentiated by the adopted measure for assessing the robustness of the schedule. In the first method [13], it is the sum of the free slack (FS) of individual activities. In the second method [14], it is the minimum value of the free slack of an activity or the minimum value of the quotient of the free slack of an activity and its execution time. In most literature sources, methods of ensuring the stability of the schedule use delay of the scheduled starting dates of particular activities, as opposed to the earliest possible dates of their commencement. These delays are defined as time buffers that are inserted before each activity in the predictive schedule. The role of each buffer is to absorb the extensions of the predecessors' duration [15]. Time buffers are assigned to individual activities, with the size being a fraction of a total slack of a given activity. The fraction of the slack is determined by looking for a compromise between the risk of failure to meet the scheduled completion date of the project and the discounted cost of its implementation. Risk assessment of delays and the discounted cost of the project's implementation are calculated using simulation techniques. An alternative to this approach is to determine time buffers on the basis of minimizing the expected cost of schedule instability [16,17]. Assuming that untimely commencement of activities generates costs for the contractor, the expected cost of schedule instability is determined as the sum of the expected costs of deviations between the starting dates of activities, forecast on the basis of the simulation results of the buffered predictive schedule, and the dates of starting activities set in this schedule. When creating the simulation scheme, it is assumed that buffered activities will be started during the project implementation no earlier than on the dates specified in the buffered predictive schedule (railway policy). Therefore, delayed commencement of activities is considered to be a delay in the date of its commencement in relation to the scheduled date, taking into account the buffers. The predictive schedule is based on the expected activities' durations and is prepared in two stages. In the first step, an initial schedule is created that provides the earliest starting dates of activities determined as a result of addressing the problem of over-allocation of renewable resources.
In the second stage, a predictive baseline schedule is created with the dates of starting activities scheduled with the buffers. The solution of the activity of minimizing the expected cost of schedule instability is obtained as follows: • A simulation of a non-buffered schedule is performed, and an activity is selected with the highest expected cost of delaying the starting date due to propagation of disturbances; • A single-unit time buffer is introduced into the schedule, delaying the scheduled start of a selected activity in relation to the earliest possible date of its commencement; • The simulation is rerun, assessing the impact of inserting a single-unit time buffer on the change in the expected cost of schedule instability in the event of disruption propagation.
If increasing the buffer for the scheduled starting date of an activity no longer reduces the expected cost of schedule instability, the next activity to be buffered is selected, or the calculation is completed. The advantage of the method is the determination of time buffers ensuring that the maximum value of the random variable T of the project completion date does not exceed the directive deadline T d . On the other hand, the disadvantage of the method is the time-consuming calculations caused by a large space of feasible solutions for this problem. To overcome this inconvenience, a priority list of buffered activities can be created using the STC heuristics-starting-time criticality [16] or heuristics CIW-cumulative instability weight [18]. The criticality of the starting date of a given activity is determined as a product of the unit cost of its untimely commencement and the probability that this activity will not commence on the scheduled date. The cumulated cost of instability of a given activity can be calculated according to the following equation: where k i is a unit cost of starting time delay for activity i, and {Succ i } is a set of activities succeeding activity i. Priority in allocating time buffers, protecting the scheduled dates of starting activities, are assigned to activities with higher criticality of starting dates (STC heuristic) or activities with a higher cumulative cost of instability (CIW heuristic). It should be noted that [19] presented a method of allocating time buffers in a robust predictive schedule using the criterion of the expected cost of schedule instability. In this case, the size of the buffer delaying the scheduled starting date of a given activity is determined analytically by allocating the total slack of the sequence of activities in the initial schedule. Simulation experiments are used to determine what part of the total slack in the sequence of activities should be allocated to a given activity in the form of a buffer delaying its scheduled starting date and to assess whether the introduced buffers reduce the expected cost of the schedule instability.
It should be noted that construction works are often executed by the general contractor in cooperation with specialized subcontractors. According to the contractual obligations of the general contractor to the client and the subcontractors, there is a need to stabilize both the starting dates of individual activities and the date of completion of the entire project. This can be achieved by the predictive scheduling method, having the appropriate mathematical instrumentation. However, it is the CC method that is promoted by PMBOK as a relatively simple tool, taking into account the psychological aspects of project management and transforming a stochastic problem of variability in activity durations into a deterministic one. Unfortunately, this method proved to be too simple and has some shortcomings in terms of buffer sizing [20,21]. There are no clear principles of grouping non-critical activities into feeding chains and sizing the feeding buffers. Moreover, if the buffer size is greater than the total slack, then its inclusion in the schedule does not make sense. It is also worth noting that most of the studies on the CCM include simple project networks with single-type feeding chains. Therefore, no interactions between multiple feeding chains are considered. Such interactions can have a large impact on the propagation of disturbances on the critical chain. The new approach presented in this article is aimed at solving these phenomena.

The New Approach
In order to use the advantages of the CCM and predictive scheduling method, the authors developed a new approach to the critical chain method. This new approach assumes grouping of non-critical activities with the same total slacks in separate feeding chains and the simultaneous analysis of both classic feeding buffers at the end of each feeding chain and buffers securing the dates for starting feeding chains (see Figure 1). The latter type of buffer is especially important if subcontracting works are involved, or the contractor shifts work teams from one construction project to another. It often happens that during construction (e.g., of a wind farm) the subcontractor must start work on the date based upon the schedule and specified in the contract. This is especially important due to the high demand for specialized engineering works such as installing wind turbines, etc. In such cases, the buffer will increase the probability that the subcontractor will start the works on the agreed date. Moreover, the use of this type of buffer prevents shifting activities in feeding chains to their earliest dates (which, as stated above, are not always desired by the planner). To distinguish the terminology, authors called these buffers:

•
Inflow buffers (IB), securing starting dates of feeding chains; • Outflow buffers (OB), i.e., classical feeding buffers securing critical chain from disturbance propagation. In order to best meet their intended purpose, inflow buffers should have a size of at least 1 time unit.
It is important for the planner/manager to decide how to manage the starting dates of activities in feeding chains. This decision is crucial for the development of the proper calculation model for the simulation of the probable scenarios of the project course. In the case of critical chains, there is no such problem since it should happen according to the ASAP rule (as soon as possible). Additionally, according to the ALAP rule (as late as possible), the first activity in a given feeding chain should start not earlier than the planned date (railway policy). However, this policy is not so obvious in relation to the other activities in a given feeding chain. For example, it is not reasonable to apply railway policy if these activities are to be carried out by the same working teams.
In addition, the aim of the authors in conducting the research was to choose the appropriate method of buffer sizing. The analytical methods used thus far, e.g., the cutand-paste method or the root-square-error method, lead to establishing relationships between buffer sizes resulting from differences in the chain's lengths and the differences between "pessimistic" and "aggressive" estimates of the activities' durations. It is not proven that, for any given network topology, having these relations between the buffers' sizes results in the highest probability of completion of a construction project on the scheduled date than if these relations are not considered. Bearing in mind that simulation is the best way to verify the effectiveness of the buffer sizing method, the essence of the approach proposed by the authors is the two-stage buffer sizing with the use of deterministic optimization and stochastic optimization techniques. Using a deterministic optimization, the initial sizes of outflow buffers are determined. Next, using stochastic optimization, these initial sizes of outflow buffers are adjusted with regard to the various probable scenarios of a project course. Determining the initial sizes of outflow buffers enables stochastic optimization to take less computing time. The optimization criteria include the construction financial or economic parameters, e.g., the monthly demand on contractor's capital for project financing CF max or the economic value of the construction project NPV, in the case of deterministic optimization, and E(CF max ) or E(NPV) in the case of stochastic optimization. The decision variables in the optimization are the sizes of outflow buffers, varying from zero up to the total slack of each activity in each feeding buffer. The sizes of the inflow buffers are calculated upon the given solution of the optimization as follows: • The size of the project buffer PB where • T-scheduled construction completion date; • T d -the completion date set by the client.
The size of the given inflow buffer IB l In reference to the suggestion of [14], the criterion of the deterministic optimization could be, for example, the maximization of the minimum of the free slacks of non-critical activities as follows: Max where • es j -the earliest start of activity j; • e f i -the earliest finish of activity i.
However, from the contractor's point of view, the economic and financial parameters of the project seem to be of greater interest. It is impossible to carry out a project with the capital requirements exceeding the financial ability of the contractor, and it is unreasonable to carry out a non-profitable project. Therefore, it is more appropriate to base the optimization criteria on the construction project cash flow analysis. According to the priorities of the contractor, the criteria of construction project cash flow optimization may include [22][23][24][25][26][27] the following: (a) Minimization of the monthly demand on contractor's capital for project financing where CF − n is the negative cash flow at the end of the current time period (month); (b) Maximization of the economic value of the construction project The solution is the buffered schedule developed in accordance with the adopted deterministic optimization criterion and meeting the following constraints: (1) The starting time of the construction project should be on a zero date (2) The starting times of activities should be non-negative (3) The relations between the starting times of activities should match the type of relations between activities; for the finish-to-start relations, the following constraint is applicable: (4) The consumption of renewable resources matches their availability ∑ p∈{A(t)} r kpt ≤ R k ; k = 1, . . . , K; t = 1, . . . , T; where • {A (t)}-the set of activities, performed on day t • r kpt -consumption of the kth resource to execute the activity p on day t; • R k -availability of the kth resource; • K-number of types of renewable resources for the execution of activity p; • T-scheduled construction completion date.
(5) The size of each outflow buffer should be a non-negative integer number where l is the number of a given feeding chain; (6) The size of each inflow buffer should be a non-zero integer number The scheduled construction completion date T should not exceed the completion date T d set by the client In the stochastic optimization stage, adjustments are made to the buffer sizes in a way to prevent the propagation of disturbances between the individual feeding chains and between the critical chain and the feeding chains. Optimization criteria include the parameters E(NPV) or E(CF max ), determined on the basis of a cash flow analysis.
A detailed description of the proposed approach is presented below.

The Proposed Procedure
The proposed approach covers the following steps of the critical chain method: 1.
Preparation of the balanced initial schedule; 2.
Preparation of the baseline schedule; 3.
Preparation of the buffered schedule.
The starting point is to model the construction project in the form of a directed, acyclic and coherent activity-on-node (AON) network, with one initial activity j = 1 denoting the start of construction and one final activity j = J denoting completion of construction. The individual activities are denoted by the condition i < j, where i is the symbol of the activity preceding the activity j. Based on the analysis of the network model with the "pessimistic" estimates of activities' durations, an initial schedule is created. According to the classical CCM, the latest starting dates of activities are taken into consideration. Then, due to the resource leveling with respect to the resource constraints, a balanced initial schedule is obtained. Resource leveling is considered here as a search for such a set of additional resource relations between the works that will deconcentrate the works to meet the planning objective (the minimization of construction project makespan) while taking into account the existing planning constraints (technological relations between works, the permissible level of consumption of individual resources, etc.). This can be achieved with the use of priority rules. As a result of using a selected priority rule, the scheduled latest date of starting work j, which has lower priority, is delayed in relation to its latest starting date, shown in the initial schedule as resulting from only technological relations between works. Such delay results from the fact that the original set of direct predecessors of activity j will be affected by a specific activity that gained priority in resource allocation over work j (in accordance with the agreed priority rule). Therefore, new resource relations are added in the network model of the project without the violation of the existing technological relations. For the identification of additional resource relations between activities, the authors proposed a method that assumes the use of binary variables x n as follows: • For x n = 1: there is a need of inserting additional resource relations regarding predeceasing and succession of works; • For x n = 0: there is a lack of such necessity.
Therefore, if N is a number of all possible additional resource relations agreed during the selection process of a priority rule, then the number of all possible solutions equals 2 N . To generate values of binary variables x n and find an optimal solution, a metaheuristic algorithm can be applied, working according to the chosen strategy of searching the N-size range of results. This algorithm solves the problem by checking if inserting new, additional resource relations satisfies the planning objective (the minimization of construction project makespan) while fulfilling the resource allocation constraints.
The authors decided to use the Tabu Search (TS) algorithm. Its advantages were proven in many scientific publications [9,10]. The TS algorithm searches the solution space by a sequence of programmed iterations. During this process, some steps are considered tabu moves, which means that they are considered to be forbidden. Due to the capability of storing the information about previously checked solutions in form of tabu lists, the algorithm avoids being stuck in local optima. Such lists are growing as the algorithm proceeds. When they reach their maximum capacity, the oldest entries of the tabu list are being erased.
In the second step of the method, the baseline schedule is obtained by reducing the activities' durations from "pessimistic" to "aggressive" values. This schedule is balanced yet not buffered. For the development of the buffered schedule, the critical chain and feeding chains have to be identified and sized. As in the original CCM, the critical chain is identified as the resource-constrained critical path in the baseline schedule. Feeding chains are identified by the grouping of non-critical activities with the same total slacks. Next, buffers are initially sized and then adjusted due to the two-stage construction project cash flow optimization. Two calculation schemes of activities' starting and finishing dates are used. The calculation scheme for the deterministic optimization is based upon the backward analysis of the project network. The calculation scheme for simulation in the stochastic optimization is based upon the forward analysis of the project network and executed in accordance with the specific decision on the policy of starting dates of activities in feeding chains. It should be noted here that if the uncertainty of the duration of a given activity is modeled by the probability distribution, usually the right tail of the probability distribution is very long (in CPM literature, pessimistic durations are twice as long as aggressive durations). However, this approach seems to be purely theoretical. An experienced construction planner is aware (to some extent) of the conditions on the specific construction site. It is the common contractual practice to cover not foreseen circumstances (defined sometimes as the acts of God) by the specific contract provisions. Therefore, similar to the statistical control of production processes, the duration of a given activity (to complete on a given construction site) may be considered by the planner as a random variable with a lognormal or symmetrical triangular distribution. Such approximation is good enough, especially when planners must base the estimation on their own experience rather than on hard-to-obtain historical/empirical data.
Finally, the buffers are inserted into the project baseline schedule. The results for an example schedule are presented in the next section.

Case Study-Wind Farm
The analyzed project is a construction of a small wind farm consisting of three wind turbines, each turbine with a 2.5 MW total power output. Each of them has a hub height of 99.5 m and a rotor diameter of 113.0 m (for a total height of 156 m). Additional major components of the analyzed project include step-up transformers located near the base of each turbine (step-up voltage: 0.69-27.6 kV), turbine laydown and storage areas, underground fiber optic cabling, and electrical collection lines (27.6 kV), with ancillary equipment, substations, and access roads. The wind turbines are spread along the terrain, resulting in slightly different conditions (including activity durations) for each of the turbines. Turbine installation activities for each wind turbine were split into two parts in accordance with the technology of hub and rotor installation. The original/initial schedule and network of the project are presented, respectively, in Figures 2 and 3.  In the case of the analyzed project, only workers employed by the general contractor (own workers) were treated as a resource and were monitored because other resources were provided by subcontractors. The maximum allocation of the resource was 54 workers (employed by the general contractor). However, the initial schedule required 56 workers (see Figure 4); therefore, the original schedule needed improvements. In the next step, the original schedule was subject to resource leveling. For the given case, the authors used a new method, as described in Section 2.1 The Method. The leveled schedule was obtained with the use of a standard, deterministic version of the Tabu Search metaheuristic algorithm. The new, leveled schedule is presented in Figure 5. The network for this schedule, including additional resource relations, is presented in Figure 6. The leveled resource graph is presented in Figure 7.   According to the proposed method, after the schedule leveling process was completed, critical and feeding chains were identified. Then, buffers were introduced to the new schedule: outflow buffers (OB), inflow buffers (IB), and project buffers (see Figure 1).
It was assumed here that the purpose of optimization was to design two variants of wind farm construction schedules, giving alternately:

•
Minimizing the general contractor's monthly demand for funds intended to cover expenses related to financing the construction investment (5); • Maximizing the economic value of the project from the point of view of the general contractor (6).
The authors used a deterministic version of the Tabu Search metaheuristic algorithm, to calculate the sizes of time buffers. Each optimization was performed both for the case with fixed relationships between buffer sizes (calculated with the RES method) and without considering such relations. A total of 5000 tests were run in order to check the correctness of the calculations. Results were compared to the RES method and are presented in Table 1. The schedule buffered with a deterministic optimization was later subject to simulations in order to compare the efficiency of deterministic and stochastic buffering optimization. Therefore, the authors included standard deviation (Std dev) and probability (P) in all result tables. Results with no relationships (new method) were better in terms of the standard deviation of makespan and probability of the project not exceeding the deadline.
The authors compared optimization results for both minCFmax and maxNPV, in terms of the classic method, with relationships between buffer sizes, and with the use of the proposed approach described in this article (Table 2). The results obtained with the use of the proposed approach for both minCFmax and maxNPV optimization show unambiguous improvement of the objective function.
In the stochastic step, solutions obtained in the deterministic model were enhanced and further developed with stochastic optimization with the use of activity durations being represented by probability distributions. The results of 5000 simulations (for each column) are presented in Table 3. The authors added additional constraints, assuming that the probability of finishing the project in contractual time is "not lower than". As presented in Table 3, minE(CF max ) optimization case delivered results with a low expected probability of finishing the project within the deadline T, which might not be satisfactory for the client. An additional constraint was introduced ensuring that all optimization results with the probability of finishing on time being lower than 70% were deemed infeasible (P(t ≤ T) ≥ P 70 ).  Figures 8 and 9 present, respectively, the best scheduling solutions obtained for max E(NPV) and minE(CF max ) objective functions.  The results obtained with the use of the proposed approach for both minCFmax and maxNPV optimization show unambiguous improvement of the objective function. The size of the difference in results is admittedly not large, but it results from the applied strict distributions for the times of performing activities. In the practice of building wind farms, the times scheduled for the performance of individual activities may differ from the actual ones usually by a maximum of 10-15%. Therefore, such distributions were defined in this case study. Meanwhile, in the literature on CCM, these differences amount to 50-100%, which is completely inconsistent with practical problems.

Discussion
The obtained results show an important improvement of the indicators in relation to the classical methods used. The proposed approach can be used, inter alia, to optimize the construction process of investments using renewable energy (such as wind farms).
The new approach uses the idea of grouping non-critical activities with the same total slacks in separate feeding chains while simultaneously analyzing both classic feeding buffers at the end of each feeding chain and buffers securing the dates for starting feeding chains (inflow buffers). It was proved that the latter (new) type of buffer is especially important. It can be used, for example, if subcontracting works are involved, or the contractor shifts work teams from one construction project to another. According to the presented research, the new inflow buffers increased the probability of successfully finishing the project. Moreover, the use of this type of buffer prevents shifting activities in feeding chains to their earliest dates (which, as stated before, is not always desired by the planner).
Interestingly, the results clearly show the NPV and CF optimization assumptions. In the case of NPV, after the buffer sizes were freed (no relationships), the buffers increased, shifting the activities to the earliest possible date. In the NPV method, it generally pays off to accomplish the required tasks and activities as early as possible.
The situation was different for the CF indicator. In this case, the scattering of construction works was beneficial, i.e., a reduction in the number of works performed in the same time period.
Another parameter analyzed by the authors was the standard deviation of obtained results. The analysis showed a great potential of the method and repeatability of obtaining valuable results.
The authors plan to develop the method further, as the results so far have been very promising and may allow solving the problems with the selection of the buffer-sizing technique. In addition, the use of the proposed approach protects non-critical pathways against disturbance propagation better than the classic critical chain method.
Author Contributions: Conceptualization, methodology, software, validation, formal analysis, investigation, resources, data curation, writing-original draft preparation, writing-review and editing, visualization, funding acquisition: J.K., N.I., J.R., and J.Z. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding.