Fractional modeling of urban growth with memory effects

According to the report provided by the United Nation’s Department of Economic and Social Affairs, the world’s population living in urban areas is expected to increase from 55% in 2018 to 68 $%$ by 2050. A better understanding of city growth will help to plan for the sustainability of cities in terms of their key attributes, such as infrastructure development, transport planning, resource allocation/consumption, and many others, which depend on the growth of population over the years. Multiple models and scaling laws have been constructed in order to understand this crucial dynamics of city growth.^1–10 The complexity of cities¹¹ is an interesting topic studying various relationships between the sizes (typically characterized by the population of the cities) and the dynamical attributes of the cities. For example, the income segregation experienced in places and by individuals can differ greatly even within close spatial proximity.¹² The complexity science of urban dynamics and its scaling laws can be referred to the two books by West¹³ and Bettencourt.¹⁴ 

Based on the model,¹ one may extract the characteristics or parameters of the cities, which can be difficult to be measured directly. For example, parameters, such as the maintenance cost to sustain the development of the city, are important for sustainable urban growth.^12,15,16 Recently, the time-dependent behaviors of urban dynamics have also been studied.¹⁷ These models, limited by the use of integer order derivatives, do not account for any memory effect. However, it is natural that the past experience of an individual living in a city influences the decision-making process. The presence of a memory effect implies that the specific system dynamics may depend on the historical states of the system more. Under the assumption that the past history of a city will influence the evolution of the urban growth at time $T$⁠, the city may respond to the history of the city on a finite interval before $T$ (memory effects) with more influence as compared to a recent and instantaneous event happening at T (no memory effects). This assumption suggests us to include memory effects with power law fading,^18–22 which can be described using fractional derivatives of non-integer order. Fractional calculus has emerged as an effective mathematical tool to model a time-evolutionary system with a memory effect in areas, such as economics^{20,21,23–29} and biology,^19,30–36 among others. Mathematically, we can incorporate the memory effect into the ODE model¹ by replacing the differential time-operator with a fractional differential operator.

In doing so, a new parameter $0<α≤1$ is introduced to capture the degree of attenuation of memory about the historical changes. Smaller values of $α<1$ imply more memory effects, and $α=1$ recovers to a non-memory model. Its exact numerical value can be determined by matching our proposed model to any given empirical attribute as a function of population over the years. If the comparison shows $α=1$ to be the best fitting, this implies that the memory effect is not important for this attribute, and it recovers the normal model.¹ If the best fit is $α≠1$⁠, it suggests the importance of a non-locality effect, which will be studied in detail in this paper. Thus, this fractional city growth model will provide additional capabilities beyond the prior models. Note that $α$ is an unique and data-driven parameter, which will be determined for different attributes associated with various cities or countries based on its given empirical scaling law between the attribute and the history of the population growth.

The remainder of the paper is organized as follows. The original urban growth model referred to as an integer order model (IOM) is introduced briefly to show its generalization to a fractional order model (FOM). Procedures on using the empirical data together with IOM and FOM are described. In our experiments of many test cases, we have explored different countries and attributes, such as congestion delay of 101 USA cities, water supply of China cities, and electricity consumption for selected countries (Singapore, New Zealand, Canada, and Switzerland). Time-dependent power scaling of $β(t)$ will also be introduced to show its improvement over some test cases as compared to constant scaling. Using the theoretical power scaling law proposed by the network interaction model,² we further study FOM with $β=5/6$ (sub-linear scaling) and $β=7/6$ (super-linear scaling) for 42 countries. Finally, the paper concludes with a summary and discussion.

To understand the behaviors of a fractional urban growth model, solutions of Eq. (5) are plotted for two different regimes: (i) sub-linear scaling with $β=5/6$ and (ii) super-linear scaling with $β=7/6$ for different values of $α$ = 0.2. 0.4, 0.6, 0.8, and 1. The figures are plotted in normalized time $t ¯$⁠, and the initial population is set to $N(0)=1$⁠. For sub-linear scaling (⁠ $β=5/6$⁠), when $N(0)< ( R / Y 0 ) 1 / ( β − 1 )$⁠, the population grows to a carrying capacity of $N ∞= ( Y 0 / R ) 1 / ( 1 − β )$ [Fig. 1(a)] at $Y 0/R=1.1$⁠. In contrast, when $N(0)> ( R / Y 0 ) 1 / ( β − 1 )$⁠, the population decreases exponentially to the same carrying capacity $N ∞$ [Fig. 1(b)] at $Y 0/R=0.9$⁠. The findings suggest that the memory effects for $α<1$ (dashed lines) will delay the population dynamics (solid line), and a longer time is needed for the population to increase (or decrease) to the carrying capacity as shown in the figures.

According to the IOM model,¹ for super-linear scaling (⁠ $β=7/6$⁠) which the growth is driven by innovation and wealth creation, the population has an exponential growth when $N(0)< ( R / Y 0 ) 1 / ( β − 1 )$ [Fig. 1(c)] and the growth will eventually reach a singularity at critical time, $t c$⁠. However, due to the limited resources, this super-linear growth is unsustainable, and it will eventually lead to a collapse when $N(0)> ( R / Y 0 ) 1 / ( β − 1 )$ [Fig. 1(d)]. If the memory effect is included for $α<1$⁠, we can see that the time taken to reach the singularity is delayed, and it will also further slow down the possible collapse of the population.

In Fig. 2, the dependence of $Y 0/R$ on the population growth is illustrated at $α=0.9$⁠. Consider a city growing with sub-linear scaling $β=5/6$ and a fixed $Y 0$⁠, a decrease in $R$ implies a reduction in its overall cost for the population growth. This allows the population to grow at a faster pace to reach its carrying capacity [Fig. 2(a)]. In contrast, for a city having its population decrease in super-linear scaling $β=7/6$ due to the limited resources, the increase in maintenance cost (⁠ $R$⁠) will speed up the depletion of resources and, therefore, accelerate the collapse of the city population [Fig. 2(b)]. The IOM cases (⁠ $α=1$⁠) are plotted (solid lines) for comparison to the specific $Y 0/R=1.8$ and 0.96, respectively.

Since $E/R$ is interpreted as the maturity time scale for an individual contributing to the urban growth, the range of $E/R$ is constrained to a reasonable range of 15–80 years, which is considered a typical working lifespan of an individual. It is also important to note that for a nonlinear optimization problem, many local minima may exist and our optimized solution is aimed to search for the global minimum. To do so, we have used a MATLAB implementation of a multi-start algorithm to help reduce the risks of being trapped in the local minima. As an extra precaution, every nonlinear optimization has been performed repeatedly for at least three times to confirm that the obtained solution is convergent and stable.

In this section, we first utilize a set of congestion delay data compiled by the Texas A&M Transportation Institute,⁴⁴ which provides the total delay (measured in person-hour) caused by traffic congestion for 101 USA based metropolitan statistical areas (MSA) from year 1982 to 2017. In Fig. 3(a), the annual congestion delay is plotted against population for all the 101 MSAs over 36 years. The overall scaling is $β$ = 1.3651 (with $R 2$ = 0.8690), which suggests a super-linear scaling (⁠ $β>1$⁠) for congestion delay because more people will lead to longer delays if the infrastructure is not upgraded. Segregating the data according to years shows yearly $β(t)$ = 1.2–1.45 from 1982 to 2017, with a variation of about 10 $%$ with respect to the overall $β$ = 1.3651.

To compare the accuracy of fractional (FOM) and non-fractional (IOM) models, the overall $β$ = 1.3651 is used to perform the matching between the empirical population $N(t)$ and the calculated results from Eq. (5). A threshold of 5 $%$ MAPE (mean average percentage error) is chosen to measure the degree of accuracy with the matching, where MAPE $≤$ 5 $%$ is considered good fitting. In Fig. 3(b), we observe that FOM (⁠ $α<1$⁠) has more smaller MAPE cases than IOM (⁠ $α=1$⁠), i.e 87 vs 73 cases. For large MAPE cases, FOM has less numbers too (14 vs 28 cases). Thus, FOM (⁠ $α<1$⁠) is found to be better than IOM (⁠ $α=1$⁠) in fitting to empirical data. A histogram of $α$ for the 87 MSA (with small MAPE) is presented in Fig. 3(c), which indicates that the best fitted $α$ is distributed uniformly in the range of $0.2<α<0.9$⁠. The corresponding values of the optimized parameters (⁠ $α$⁠, $R$⁠, and $E$⁠) for each MSA with good fitting MAPE $≤5%$ are shown in Tables S1 and S2 of the supplementary material.

In Figs. 3(d)–3(f), we show three selected MSA examples with different MAPE values: $≪5%$⁠, $≈5%$⁠, and $≫5%$ MAPE to provide a visual correspondence of the degree of fitting. In Fig. 3(d), we observe an excellent agreement between both models and population growth of Worcester, with FOM (⁠ $α=0.7741$⁠) having smaller MAPE (0.6 $%$⁠), which is better than IOM (1.4 $%$⁠). The estimated maintenance cost (⁠ $R$⁠) and the additional cost (⁠ $E$⁠) using IOM and FOM are summarized in Table I. For this case, the estimation of $R$ and $E$ using IOM and FOM is relatively close, with FOM estimating a 10 $%$ more maintenance and additional cost compared to IOM.

For cases with average MAPE around 5 $%$ [see Fig. 3(e), MSA is Charlotte], noticeable discrepancies are observed, but FOM (⁠ $α=0.5258$⁠) remains better than IOM (MAPE = 4.8 $%$ vs 6.6 $%$⁠). The estimated maintenance cost (⁠ $R$⁠) and the additional cost (⁠ $E$⁠) using IOM and FOM are summarized in Table I. The estimated $R$ and $E$ using FOM are higher than IOM. Considering that the MAPE of FOM is lower and in good agreement with the population data, IOM may have underestimated the maintenance and additional cost.

For cases with large MAPE $≫5%$ [Fig. 3(f)], IOM and FOM are unable to provide good fitting. After filtering out cases (from the total 101 cases) with MAPE larger than 5 $%$⁠, we have 74 and 87 cases, respectively, for IOM and FOM, which confirms again that FOM is better than IOM in providing better matching. In each of these 87 cases, MAPE of FOM is always smaller than that of IOM.

In this section, we examine the urban data of China from year 2004 to 2019 available from the China Urban Construction Statistical Yearbook. The data contain various urban attributes of over 600 cities (including prefecture-level and county-level) in China. Among the attributes, water supply, which is highly relevant to urban sustainability, is chosen to be studied in more detail. Power-law fitting is performed to extract the scaling exponents $β$ of water supply with respect to the population at the country level as well as the province level. Figure 4 shows $β=1.1274$ (country level) and $β=1.0047$ (Hebei province level) for water supply. Using these two values of $β$⁠, both IOM and FOM are applied to analyze the urban population data for Hebei province. There are 28 cities (some cities are excluded due to lack of data) in Hebei province, of which the population growth is modeled with $β=1.1274$ (country level) and $β=1.0047$ (province level) extracted from the water supply data (see Fig. 4).

A comparison of the calculated MAPE for IOM and FOM is shown in Figs. 5(a) and 5(b). Out of the 28 cities, nearly all MAPEs from FOM are equal to or less than IOM except for one city (Xingtai) and thus confirms the improvement of FOM over IOM. There are ten cities that FOM recovers to the IOM limit showing that $α=1$ is the best fitting. The corresponding $α$ for the other 17 cases is in the range of $α=0.2$–0.7, plotted as black circles in the figures. It is important to note that for some cities, such as Xinle, Huanghua, and Zunhua, the improvement of FOM is significant—from large MAPE of more than 10 $%$ to less than 5 $%$⁠. The optimized values of parameters (⁠ $α$⁠, $R$⁠, $E$⁠) for both IOM and FOM with $β=1.1274$ and $β=1.0047$ can be found, respectively, in Tables S3 and S4 of the supplementary material. Using a threshold of MAPE = 5 $%$ to assess the degree of good fitting, the adoption of FOM has confirmed its advantages over IOM. As examples, Nangong and Zunhua from Hebei province are presented in Fig. 6 for comparison: solid lines (FOM) and dashed lines (IOM).

The details of the improvement of FOM over IOM can be referred to the figure caption. The estimated maintenance cost and the additional cost for Zunhua based on IOM and FOM using different value of $β$ are summarized in Table II. The improvement is more significant for a large MAPE case obtained from IOM as shown in Fig. 6(b).

In this section, we will test the models in using the electricity consumption data of four countries: Singapore, New Zealand, Canada, and Switzerland. The electricity consumption data of Singapore between year 2005 and 2018 are available from the Department of Statistics. The data for the other three countries between year 1990 and 2019 are obtained from the International Energy Agency. A comparison of empirical data fitting using FOM and IOM on population growth of the four countries is shown in Fig. 7 with their $β$ values plotted in the insets.

In general, FOM again confirms that it has improved over IOM in having better agreement with smaller MAPE as shown in the cases of Singapore and New Zealand. Contrast to IOM, FOM is able to capture the population growth of Singapore and New Zealand more accurately, further reduces the MAPE by nearly half from an already small MAPE of 3.4 $%$ (IOM for Singapore) and 1.57 $%$ (IOM for New Zealand). In the cases of Canada and Switzerland, FOM does not exhibit additional advantages over IOM and thus recovers to the IOM. The details of the calculations, such as MAPE, $R$⁠, $E$⁠, and $α$⁠, between IOM and FOM can be found in Table III.

While it is difficult, if not impossible, to determine $E$ and $R$ in practice, we could reference them to the average electricity consumption per capita for estimation of the order of magnitude. We estimate the electricity consumption per capita based on the collected data and plotted in Figs. S1–S4 of the supplementary material. It is evident that the cost for electricity consumption evolves as a function of time. Here, it is important to emphasize that the values of $E$ and $R$ calculated are the characteristic values over the years, i.e., 2005–2018 for Singapore and 1990–2019 for the other 3 countries. Based on the collected data, the average electricity consumption per capita is determined to be 8.40, 9.17, 16.53, and 7.77 for Singapore, New Zealand, Canada, and Switzerland, respectively. The average costs are of similar order of magnitude to the estimated $R$ by IOM and FOM shown in Table III. Considering that FOM provides better matching to the data in the cases of Singapore and New Zealand, IOM may have underestimated (for Singapore) and overestimated (for New Zealand) the maintenance cost, $R$⁠. For Canada and Switzerland, FOM essentially recovers to IOM, and both models provide the same estimation for the maintenance cost, $R$⁠, and additional cost, $E$⁠.

To explore the improvement of using time-dependent $β(t)$ scaling for both FOM and IOM, we apply the testings to the traffic congestion delay of USA cities studied above. Apart from obtaining $β(t)$ in an average sense over the entire dataset, noted as $β 0(t)$⁠, it is also possible to segment the attribute dataset by cities and further extract information more specific to an individual city, noted as $β i(t)$⁠. This allow us to study if such localized information is helpful for the modeling. We decided to investigate the two types of time-dependent $β(t)$⁠: $β 0$ and $β i$⁠. For first type, data fitting of Eq. (8) is performed on the entire set of 101 MSAs and the resulting scaling is denoted as $β 0$⁠, which is the $overall$ time-dependent $β(t)$ extracted from data containing ALL cities. In doing so, the $R 2$ has improved slightly from 0.8690 to 0.8800. For the second type, the data fitting of Eq. (8) is performed ONLY on individual city and the resulting scaling is denoted as $β i$⁠, which is the $individual$ city-based time-dependent $β(t)$⁠. In doing so, the $R 2$ for the individual city is larger than 0.85 with the majority in the range of 0.95–1.0. In comparison, the original constant scaling is denoted as $β c$ in the subsequent discussions.

IOM and FOM for both $β 0$ and $β i$ are applied to model the population growth of 101 MSAs. With the same criteria of MAPE $=5%$ chosen to assess the degree of good agreement, Fig. 8(a) shows a comparison of the number of good agreements (⁠ $MAPE≤5%$⁠) for different types of models (IOM or FOM) in using different time-dependent scalings, $β 0(t)$⁠, $β i(t)$⁠, or constant $β c$⁠: (a) IOM and FOM with $β c$ and (b) t-IOM and t-FOM are IOM and FOM with different time-dependent $β(t)$⁠. It is observed for the same type of $β$ used, and FOM is always better than IOM. Among all three models of FOM, t-FOM using the overall time-dependent $β 0$ has the high number of good agreement (96 cities) as compared to FOM using constant scaling (87 cities). Interestingly, t-FOM using individual time-dependent $β i$ is significantly under-performed with only 46 good agreements. Similar findings are observed for IOM cases, and this suggests that the overall time-dependent scaling $β 0(t)$ is probably the better approach. The combination of FOM using the overall time-dependent scaling $β 0(t)$ has provided the best approach in having 96 (out of the 101 MSAs) cases with better good agreement of $MAPE≤5%$⁠. As shown in Fig. 8(b), with the incorporation of time-dependent scaling of $β 0(t)$⁠, the estimated values of $α$ have shifted toward the higher values with majority in the range of 0.9–1.0.

In Figs. 9(a) and 9(b), two MSAs: Phoenix–Mesa and Nashville–Davidson are selected to illustrate the modeling of population growth in using FOM with their time-dependent $β 0(t)$ and $β i(t)$ plotted in Fig. 9(d). Phoenix–Mesa is an example of degrading fitting when $β i$ is adopted in comparison with $β 0$⁠, with MAPE = 1.49%–12.75. Nashville–Davidson is an example of improved fitting when $β i$ is adopted in comparison with $β 0$ from MAPE = 5.08 $%$ to 2.36 $%$⁠. The estimated maintenance cost and additional cost based on FOM using different forms of $β$ for Phoenix–Mesa and Nashville–Davidson are summarized in Table IV. The corresponding values of the optimized parameters (⁠ $α$⁠, $R$⁠, and $E$⁠) using the overall time-dependent $β 0(t)$ for MSAs with relatively good fitting (MAPE $≤5%$⁠) are shown in Tables S5–S7 of the supplementary material.

To probe further, we may examine the fitting of the congestion delay using Eq. (8) with $β 0$ and $β i$ on both MSAs. Compared to $β 0$⁠, it is obvious that $β i$ has improved the empirical fitting on data (circles) significantly as shown in Fig. 9(c): (I) Phoenix–Mesa: $R 2$ is improved from 0.7414 to 0.9972 and (II) Nashville–Davidson: $R 2$ is improved from 0.5337 to 0.9953. The time variation of extracted $β i$ of both MSAs is also plotted in Fig. 9(d), which shows radically different time-varying trends. For Phoenix–Mesa, we have an increasing $β i$⁠, whereas $β i$ of Greensboro is decreasing over time. The comparison shows that $β i$ of Nashville–Davidson behaves consistently with yearly $β$ and $β 0$⁠, which suggests that the time-dependent scaling of $β i$ can improve Nashville–Davidson, but not for the Phoenix–Mesa case, which shows poor agreement in using $β i$ as discussed above. In general, the overall time-dependent scaling (⁠ $β 0$⁠) has the best performance in our study with the highest number of good agreements. The individual time-dependent scaling $β i$ is better only in some selected cases, where 75 out of the 101 MSA has shown higher MAPE (poor agreements) in using $β i$⁠.

Instead of extracting the required $β$ for specific attributes, we consider the two theoretical scalings, which have been derived in Ref. 2: $β=5/6$ (sub-linear) and $β=7/6$ (super-linear). These theoretical values of $β$ allow us to apply FOM to analyze many countries independent of the specific attributes. It is expected if certain attributes of the studied countries are within the values of $β=5/6$ and $β=7/6$⁠, and the calculated results to be shown below will be valid. A random selection of 40 countries from the World Bank between year 1960–2018 was used for their population data. With additional two more countries (Singapore and Hong Kong) in using their government reported data, we have a total of 42 countries to be studied.⁴⁵ 

The calculated results are presented in four different combinations: IOM with $β=5/6$ and $β=7/6$ and FOM with $β=5/6$ and $β=7/6$⁠. Figure 10 shows the histrograms of the calculated parameters: (a) MAPE, (b) $Y 0/R$⁠, (c) Age (⁠ $E/R$⁠), and (d) $α$ calculated from FOM. A threshold of MAPE = 5 $%$ or lower is used as an indicator of good agreement. Figure 10(a) shows that majority of the countries (at least 70 percents) have good agreement in using FOM for both values of $β$⁠. For IOM in using $β=7/6$⁠, it shows worse performance with only 14 countries (one-third) shows good agreement. In all countries, FOM is always better or equal than IOM independent of $β$⁠.

From the figure, we can draw some key findings. First, we notice that the maturity age (⁠ $E/R$⁠) lies near the boundary of the range of age (15–80 years) that we used in our models assuming that these are periods of an individual to have strong effects on the dynamics of the cities. For $β=5/6$⁠, $E/R$ is mostly at $E/R=15$–20 years for both IOM and FOM. For FOM with $β=7/6$⁠, it is mostly at $E/R$ = 70–80 years. Second, the values of $Y 0/R$ are mostly within 1–2 for both IOM and FOM. Finally, the value of $α$ from FOM is mostly is about 0.8–1 (for $β=5/6$⁠) and is about 0.4–0.6 (for $β=7/6$⁠). These findings represent the general characteristics of FOM models for sublinear and superlinear scalings over a wide range of countries. The optimized values of the variables (⁠ $Y 0/R$⁠, $E/R$⁠, $α$⁠) for both IOM and FOM with $β=5/6$ and $β=7/6$ can be found in Tables S10 and S11 of the supplementary material.

Two countries are selected to illustrate the memory effect at different $α$⁠. The population growth of South Korea is plotted in Fig. 11(a) with the best fitted FOM for $β$ = 5/6 having MAPE of 1.47 $%$ at $α 0=0.6816$⁠. We use a 20 $%$ deviation from this value of $α 0$⁠, which are denoted as $α 0↑$ and $α 0↓$⁠, to study the effects of varying $α$ in the FOM results. In Fig. 11(a), it is observed that smaller $α 0↓$ will lead to faster growth than $α 0$ in the initial growth period (1960–1980) and slower growth than $α 0$ at the later stage of the growth (after 1980). On the other hand, larger $α 0↑$ will lead to opposite trends with slower growth than $α 0$ in the initial growth period but faster growth than $α 0$ at the later growth period.

The population growth of United States is plotted in Fig. 11(b) with the best fit FOM for $β$ = 7/6 having an MAPE of 0.65 $%$ at $α 0$ = 0.7906. Here, larger $α 0↑$ will always lead to slower growth than $α 0$ over the entire period, and smaller $α 0↓$ will always lead to faster growth than $α 0$⁠. The two different dynamics of memory effect at changing $α$ observed is due to the range of maturity age (⁠ $E/R$⁠) obtained in the optimization process. If we have a smaller $E/R$ near the lower limit (15–20 years), we have the memory effect of the typical South Korea scenario as shown in Fig. 11(a), which is typical. If the obtained $E/R$ is large (70–80 years), we will have the USA-like memory effect as shown in Fig. 11(b).

Limited by the ordinary calculus, a simple ordinary differential equation (ODE) model is incapable of capturing memory effects. The integer order model of Eq. (2) is applicable only on the condition that the city forgets all past experiences. This no-memory (integer order) approach cannot always be used in the socio-economic context. It is natural that an economic agent reacts based on historical events with the recent event being more influential. This work assumes that today’s population growth may remember previous notable changes in population growth over time $T$⁠. The memory effect implies that at repeated similar changes, the population growth rate can be affected differently than the normal ODE model. This leads to the fact that the same value of growth may correspond to a number of different prior time intervals. To capture this, we have applied a fractional calculus to incorporate a memory effect with power law fading into the urban growth model pioneered by Bettencourt,¹ and we have shown that this approach provides better fitting to empirical data. Thus, in our approach, the ordinary differential operator in the ODE is replaced by the Caputo fractional derivative that we term as a fractional order model (FOM).

It is observed that FOM may improve the accuracy over the normal ODE model by analyzing a large number of test cases of different cities, countries, and attributes: (a) Congestion delay of 101 cities in the USA; (b) water supply for cities in China; (c) electricity consumption of four countries (Singapore, New Zealand, Canada, and Switzerland); and (d) 42 countries in using theoretical values of $β=5/6$ (sub-linear) and $β=7/6$ (super-linear). The constant power scaling component may be improved by using suitable time-dependent power scaling law, which is empirically determined from a given set of data. We believe that this model is useful to understanding the urban dynamics with scaling law in the range of $β=6/7$– $β=7/6$⁠.

The fractional parameter, $α$⁠, which describes the degree of attenuation of the memory about the historical changes, can be different for different cities. For the same city, it can be different for different attributes of the city. Different values of $α$ actually show how strongly their growth is affected by their past. A city with $α=1$ is likely to follow Eq. (2) in future, and it would be robust against most external factors as it was in the past; therefore, a prediction is straightforward. For cities with $0<α<1$⁠, any future major happening, the growth is more likely to be affected in a less-predicted manner, hence following FOM. It could be interesting for the scholars in social sciences to identify what key events have potentially led some cities to acquire $α<1$ in some of their attributes. The model can be further tested with other data, and the extracted $α$ might be used to analyze and shed light into the resilience of different urban systems. Throughout this work, $α$ is limited to between 0 and 1 to avoid high order dynamics that is absent from the urban growth model by Bettencourt so that we can have a consistent and fair comparison with the current ODE model. Numerically, it is possible to relax this constraint in the optimization procedure so as to explore the existence of higher order dynamics in future studies and compare with empirical data.

On the other hand, analyzing the data using different mathematical approaches will be useful to enhance our understanding and to provide better prediction of the complex dynamics of the system. It could also be interesting to compare the fractional model with other techniques that incorporate time delays or short-run memory, such as delay differential equations.^46,47 Time delay is a common phenomenon that widely exists in economic systems and others. Furthermore, enriching the fractional model with time delays could be of interest to capture the effects of both long-run and short-run memory in the city growth dynamics. This combined approach in a fractional delay differential equation (FDDE) has been proposed in the literature.^48,49

Overall, we think that for systems that involve human actors, ignoring memory effect and assuming instantaneous response of a system is perhaps oversimplified. The decision-making process of humans often involves referencing to past memories or experiences. For this purpose, fractional calculus has emerged as an effective tool to capture memory effects and demonstrated potentials in many economic applications.^25,50–52

See the supplementary material for the optimized parameters ( α, R, and E) of the least square fitting performed on the data sets considered in Sec. III.

This work was supported by the SUTD SGP grant: Data Driven Design Solutions for Cities (No. SGPCTRS1804) and the USA ONRG grant (No. N62909-19-1-2047).

The authors have no conflicts to disclose.

Chun Yun Kee: Data curation (equal); Formal analysis (equal); Methodology (equal); Visualization (equal); Writing – original draft (lead); Writing – review & editing (equal). Cherq Chua: Data curation (equal); Formal analysis (equal); Visualization (equal); Writing – review & editing (equal). Muhammad Zubair: Writing – review & editing (equal). L. K. Ang: Conceptualization (lead); Funding acquisition (lead); Project administration (lead); Supervision (lead); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.