This modeling study was initiated based on the primary observations from the PR rule [5]. This was also influenced by the predicting potential for upper limb recovery 1, 2 (PREP 1, 2) algorithms [66, 67], and by a study [62]. Functional recovery post-stroke can be dichotomously differentiated based on the connectivity of the ipsilesional CST assessed within two weeks after a stroke. Reliance on alternative tracts increases if neural networks along the CST become too sparse or disconnected to convey motor commands to appropriate motoneuron pools. We aimed to describe the selective and competitive relationship between two hypothetically primary descending pathways after stroke, employing Hebb’s model, commonly used to model use-dependent neural plasticity [71].

Our model successfully describes a battery of phenomena we commonly observe in individuals affected by stroke. Above all, the model replicates that the degree of the initial anatomical damage to the CST substantially determines the direction of neural recovery (behavioral restitution or compensation). In our model, some cells (i.e. ipsilesional M1 cells) connected to the CST are destroyed and even disabled by stroke (their firing rates are 0) emulating biological systems [16, 63]. Suppose these cells can no longer function and the number of such cells exceeds a certain threshold. In that case, the neuromotor system begins to use relatively more optimized (less damaged) neural pathways to compensate for the role of those cells. The model assumes that corticospinal pathways connected to those cells are accordingly disabled by stroke. This aligns with the PR rule and PREP algorithms, which suggest that the initial anatomical capacity of the CST, measurable with TMS or imaging, determines the degree of potential recovery categorized as “fitters” versus “non-fitters” (PR rule), or “completed, notable” versus “limited, none” (PREP algorithm). If the cells affected by stroke retain some rudimentary function (i.e., their firing rates are not 0), meaning they are not entirely dead, the neuromotor system continues to regrow and strengthen them for motor execution [52], and those cells and their corresponding circuits become re-optimized. This can be described as the Hebbian-type refinement of neural circuits [51]. This cascade is analogous to the typical occurrences post-stroke. Indeed, stroke deprives some neurons of their normal metabolic substrates, causing them to lose function within seconds [53]. Those neurons are unable to maintain their normal transmembrane ionic gradients, eventually leading to apoptotic and necrotic cell deaths [30, 52]. Surviving neurons in the peri-infarct cortical areas with sufficient blood perfusion undergo active structural and functional remodeling [51]. Dendritic growth and axonal sprouting occur to rewire damaged connections or form new connections [52]. Accordingly, motor function recovers and improves; behavioral restitution can be achieved. This suggests that the spontaneous recovery mechanism and use-dependent neural plasticity revive ipsilesional cells or pathways if they survive [52].

Studies reported that individuals with milder impairment show quicker recovery than those with more severe motor impairment [11]. Meanwhile, motor improvement in individuals who are assumed to increase their reliance on alternative tracts continues throughout the first year after stroke. This suggests that it takes longer for individuals with greater impairment to adapt to compensatory mechanisms for the damaged CST and achieve functional fractionation of the alternative neural pathways, in comparison to the use of the recovered CST and ipsilesional hemisphere in individuals with mild impairment. When CST connectivity is deficient, compensatory strategies emerge, including a shift in interhemispheric lateralization towards the contralesional hemisphere and a shift in representational maps around the infarcted zone [52]. Although the neural circuit used as a substitute is relatively intact or less damaged, time is required for adaptation through activity- or use-dependent plasticity. The amount of fractionation correlates with the time since stroke onset [8]. In the model simulation, we described this phenomenon by assuming that the initial RST accessibility is very low (the firing rates of the RST cells were set near 0 ~ 0.1). The functional capability of the RST begins near 0, regardless of the degree of CST damage caused by stroke (refer to Figs. 1, 2 and 3), so it takes more time for the RST to reach its maximum functional capacity. As a result, motor improvement in individuals for whom the RST acts as the dominant descending pathway tends to be delayed in comparison to those for whom the CST is the dominant descending pathway.

The time evolutions of FMA subscores suggest that the capability of conducting the type of out-of-synergy movements slowly increases in stroke survivors with severe impairment throughout the subacute phase and into the chronic phase (refer to Fig. 4). This implies that the functional capabilities of the CST and RST increase in parallel post-stroke. This observation agrees with imaging studies showing that all tracts tend to recover simultaneously after damage [31]. Branches of the CST to motoneurons controlling the fingers can still function in individuals with severe impairment; those individuals merely show a degraded ability to extend distal joints due to weak CST connectivity or overwhelming abnormal synergies mediated by the RST [1, 35]. However, at the same time, there is the priority in optimization between the tracts. All those results imply that the relatively less damaged RST network is optimized with priority, and the functional capability of the CST may continue to improve in compensation for the residual functional capability, which the RST cannot fulfill solely.

Constraint-induced movement therapy (CIMT) is an established rehabilitation approach to promote the recovery of the ipsilesional hemisphere while suppressing contralesional motor drive post-stroke [40]. This therapy effectively reduces maladaptive neural reorganization by interfering with cortical motor expansion and reducing reliance on the less affected limb [61]. However, some have argued that CIMT can cause maladaptation through functional reinforcement of compensatory motor activity [37]. True behavioral restitution, a return to normal motor patterns with the affected limb [41], requires the recruitment and restoration of the residual ipsilesional hemisphere/CST. In contrast, training of movement within synergy patterns with the affected limb can selectively enhance the structural and functional capabilities of the RST. A study with intact primates showed that strength training with pulling movements enhanced the RST, not the CST [24]. In humans, the RST is predominantly used to excite the biceps on the paretic side post-stroke [62]. Also, strength training with grasping promotes the RST [45].

Our model did not differentiate between the types of movements that cause neural refinement or optimization of a particular tract. We assumed that activity- or use-dependent neural plasticity in our model is achieved with repetitive movements of muscles that can be excited via either the CST or RST. Our primary interest was to examine how the CST and RST interact and evolve, per the hypothesized principle “A tract that is more strongly optimized becomes optimized with priority”. Those muscles include shoulder abductors/flexors [21, 27], elbow flexors [62] and finger flexors [2]. Either the corticospinal network or the reticulospinal network is optimized with priority through use-dependent processes by repeating movements with those muscles or imagery movements, depending onthe initial CST connectivity. If the reticulospinal network is optimized with priority, it ultimately facilitates RST upregulation and abnormal synergies.

Abnormal co-activation across muscles may originate from the physical feature of the RST that branches into multiple motoneuron pools across the upper extremity and activates them together when a central command descends [4, 13, 14, 29, 65]. Even in individuals who predominantly use the RST, selective muscle activation may need to require the CST to convey motor commands, instead of the RST. Naturally, its repetition could promote the structural and functional capabilities of the CST and lead to its dominant use. Upper-limb movement away from stereotypical abnormal synergies may be a promising way to achieve behavioral restitution. Training stroke survivors to simultaneously coordinate and synchronize multiple fingers, alleviating the flexion synergy, is effective in improving finger individuation and hand dexterity [48]. This leads to a significant improvement in the FMA score, suggesting that the impairment of the upper extremity is alleviated.

In our model, the functional capability of the CST is limited by the number of cells on the ipsilesional hemisphere that survive after stroke. However, we believe that cortical representations can be reconstituted in the peri-infarct tissue through repetitions of out-of-synergy movement, like in CIMT [61]. It would be worth investigating the effect of inducing the corticospinal network to be optimized with priority in individuals with weak or no initial CST connectivity after stroke through conducting out-of-synergy movements and blocking activities that enhance the reticulospinal network.

We summarize the following four assumptions imposed on the model, which might be limitations of this study:

1) The improvement in motor function requires those in the functional capabilities of the CST and RST (the additive nature of the functional capabilities of the hypothetically primary descending motor tracts).

2) Specificity of the CST to the improvement in scoring in FMA out-of-synergy test items in both fitters and non-fitters.

3) Specificity of the CST to the improvement in scoring in FMA in-synergy test items in fitters.

4) Specificity of the RST to the improvement in scoring in FMA in-synergy test items in non-fitters.

The proposed model is based on Hebbian learning (behavioral improvement) occurring within these predefined assumptions with differential learning rates across (simulated) subjects. The model primarily relies on the thresholds on the CST beyond which the muscles can be innervated. Although we artificially divided the Minor and the Substantial CST Damage cases with a boundary of 20% CST dead cells, it does not influence our results and conclusion. If the functional capability of the CST does not increase beyond the thresholds, the neuromotor system improves the RST to innervate the muscles. For non-fitters, the functional capability of the CST does not reach the thresholds needed to execute a planned movement by innervating the corresponding muscles. Instead, the RST, which is less damaged than the CST, is enhanced to compensate for the role of the CST. They fail to achieve 80% motor recovery. Of course, we assume that all cells and connections have identical capabilities, which may be unrealistic. It is difficult to interpret that the damage size of the corticospinal network can always differentiate non-fitters versus fitters in practice.

This modeling study relies solely on pieces of empirical evidence, not directly measured data. By reproducing the time evolutions of FMA subscores, we tried to infer the motor tract in dominant use and interactive dynamics between the hypothetically primary motor tracts. We assumed that the out-of-synergy movements instructed in the FMA are mediated predominantly via the CST. The type of out-of-synergy movement may be mediated even via the RST, though its execution is weak. We cannot rule out the possibility of functional fragmentation (remodeling) of the alternative neural substrates (contralesional cortices and RST). Several studies demonstrated that structural reorganization of the contralesional cortices occurs following stroke and contributes to motor improvement, possibly promoting joint individuation [3, 6, 57]. However, studies evidenced that the functional upregulation of the RST leads to substantially reduced hand/arm dexterity (i.e. the capability of performing out-of-synergy movements) [10, 50]. Thus, our assumption is considered reasonable. Also, we assumed that different motor tracts are involved in performing in-synergy and out-of-synergy movements in non-fitters. Our previous study with large databases of FMA subscores revealed that significant asymmetry in scoring between in-synergy and out-of-synergy FMA test items was observed in individuals with severe-to-moderate impairment, whereas symmetry was observed in individuals with mild impairment [33, 34]. This study showed similar results (Fig. 4). Such a reverse relationship in capability between in-synergy and out-of-synergy movements, evidenced in FMA subscore distributions and imaging studies [10, 32], could justify our assumption of the involvements of different motor tracts in performing in-synergy and out-of-synergy movements in non-fitters.

Our computational model of tract dynamics involves several parameters: firing rate, weight and learning gain. We set the initial values of the firing rate and weight to identical random values within a specific range. The range was set as a function of the number of the alive CST cells (i.e. [\(\:\sqrt/10-0.1\), \(\:\sqrt/10+0.1\)]). However, setting a specific range does not influence the competitive and parallel relationship between the two tracts. Figure 6 shows the time evolutions of the simulated functional capabilities of the CST and RST, with the initial values of the firing rate and weight either dependent or independent of the number of alive CST cells.

Time evolutions of the simulated functional capabilities of the CST and RST, with (a) the initial values of the firing rate and weight dependent on the number of alive CST cells and the learning gains varying across subjects, (b) the initial values independent (randomly chosen between 0.1 and 0.9) and the learning gains varying, (c) the initial values dependent and the learning gains fixed, and (d) the initial values independent and the learning gains fixed

For simplicity, we set the values of the learning gains of the firing rate and weight for the CST and RST to identical random values within a particular range, respectively. The same range was applied to fitters and non-fitters. Figure 6 shows the time evolutions of the simulated functional capabilities of the CST and RST with the learning gains varying across subjects, in comparison to those with the learning gains fixed across subjects. Variability in the learning gains does not influence the competitive and parallel relationship between the two tracts. Although we set the same ranges for the learning gains in fitters and non-fitters for simplicity, the learning gains can reflect the difference in individuals’ capability for neural recovery and plasticity. The majority of individuals with strong initial CST connectivity (i.e. fitters) generally recover an average of 70–80% within 1 or 2 months after stroke, regardless of sex, age, race, and the amount of behavioral effort [46, 72, 74]. However, individuals with weak or no initial CST connectivity (i.e. non-fitters) show relatively large variability in the recovery rate [41, 70]. This large variability possibly originates from factors influencing the rate of neural plasticity, including sex, age, lesion size, fractional anisotropy asymmetry, necrotic tissue, edema, and inflammation, all of which are known to affect post-stroke recovery [22, 69]. Therefore, it is more reasonable to set wider ranges for the learning gains in non-fitters.

In the reproduction of FMA subscores, we adopted several assumptions: (1) each instructed movement of the FMA can be conducted once the functional capability of the CST or RST reaches its corresponding threshold; (2) conducting a movement fully (Score “2”) requires a greater functional capability of the dominant tract than is needed for conducting a movement partially (Score “1”); (3) the same thresholds for out-of-synergy movements are applied to fitters and non-fitters; and (4) different thresholds for in-synergy movements are applied to fitters (on the CST) and non-fitters (on the RST). In the current study, the threshold for each instructed movement of the FMA was not determined based on direct empirical evidence. We tuned the values of those thresholds to closely replicate the empirical results shown in Fig. 4, as long as those assumptions remain valid. While the values of the thresholds can slightly alter the amplitude, standard deviation, and shape of the resulting curves, the observation remains unchanged. We confirmed that no other scenarios with different threshold values are possible. Since our empirical data are from a small sample size, which may not represent the entire population, we did not focus on precisely matching the simulation results to the empirical ones.