Although there are various models of epidemic diseases, there are a few individual-based models that can guide susceptible individuals on how they should behave in a pandemic without its appropriate treatment. Such a model would be ideal for the current coronavirus disease 2019 (COVID-19) pandemic. Thus, here, we propose a topological model of an epidemic disease, which can take into account various types of interventions through a time-dependent contact network. Based on this model, we show that there is a maximum allowed number of persons one can see each day for each person so that we can suppress the epidemic spread. Reducing the number of persons to see for the hub persons is a key countermeasure for the current COVID-19 pandemic.

An epidemic disease is an area where a mathematical model has a potential power to derive how to suppress spreading of the disease. Although there are various types of models such as the famous Kermack and McKendrick model^{1} to agent-based models^{2} and age-structured models,^{3} there are a few individual-based models that contain a limited number of assumptions and can be analyzed analytically to suggest to each individual how they should appropriately behave. Thus, here, we introduce a simple topological model of epidemic disease spreading. The key ingredient is a contact network among people for each day. This model can allow the time-varying nature for the underlying network topology due to possible countermeasures enforced by each government or each individual including susceptible persons. In addition, this model includes an infection history of an epidemic disease so that the disease states can match with the current coronavirus disease 2019 (COVID-19) epidemic more closely. The theoretical arguments based on the model directly show that hub persons should be the most paramount point to consider the countermeasures,^{4} especially in the current COVID-19 condition.

COVID-19 has been spreading worldwide and producing serious problems. The number of infected persons is 16 950 407, and the number of deaths is 664 961 as on 30 July 2020.^{5} However, because there are no effective existing vaccines or medicines specific for COVID-19 until now, currently, we can maintain social distancing from people,^{6} wear a mask, and/or wash our hands. The second wave of COVID-19 is spreading in some countries, while the possibility and its scenarios for the second wave have been examined theoretically.^{7,8} In addition, more waves are expected^{9,10} because infections from pre-symptomatic and asymptomatic infectious people are sufficient to sustain the current COVID-19 pandemic.^{11} Here, our question is how many contacts per person for each day are few enough for suppressing the further spreading of COVID-19. We formulate a topological model for an epidemic disease and have a solid answer on the maximum number of contacts per day for each person.

The mainstream for the mathematical modeling of epidemic diseases is via differential equations.^{12,13} However, for example, the most popular Kermack and McKendrick model [Susceptible–Infected–Recovered (SIR) model]^{1} and the susceptible–infected–susceptible model^{12} are too macroscopic to derive countermeasures interpretable for each susceptible person to follow. To derive and/or evaluate countermeasures microscopically that we should enforce for each person’s level, we need to use an individual-based model such as agent-based models.^{2} However, such a model has not considered the temporal change of the infecting probabilities depending on the time after the infection. In addition, such a model so far has not reached a criterion that each individual person can understand and follow for suppressing an epidemic spread, which is an important point for the current COVID-19 pandemic.

On the evening of 27 April 2020, we found that Professor Michael Small posted a preprint^{14} and declared that the topology of a contact network is sufficient to understand an epidemic and take countermeasures. Instead of considering the reproduction number,^{15} which is hard to estimate from a short time series,^{16} their arguments are based on the topology of contact networks. It has been shown that the topology of the contact network is related to an epidemic spread on the network: for example, in a scale-free network, an epidemic is more likely to spread than in a regular lattice or a small-world network.^{17} A more mathematical work^{18} discusses how an epidemic spread will die out depending on the comparison between the speed of curing and the network’s spectral radius. Therefore, we were inspired by the work of Professor Small for constructing a topological model of an epidemic disease. In particular, rather than using the reproduction number, we would like to characterize the underlying time-dependent dynamics of an epidemic spread by a time-varying contact network or a temporal network.^{19–21} Due to the technical advancements, we have been able to collect time-dependent contact information, and such data have been widely analyzed.^{22–25} Hence, our model will have a good affinity with such time-dependent contact information. Thus, we intend to construct an individual-based model that (i) can describe the current COVID-19 situation more accurately by using the infection history as well as contact networks and (ii) the model can be analyzed analytically so that we can obtain a topological criterion that can be implemented by each *susceptible individual* for suppressing the epidemic spread. This topological criterion is why we call our model a “topological epidemic model.”

Our model considers an individual person $i$ in some area. Among them, we consider a time-varying contact matrix $A(t)$ at day $t$. In particular, if persons $i$ and $j$ have at least a contact on day $t$, we have $ai,j(t)=aj,i(t)=1$. Otherwise, we set $ai,j(t)=aj,i(t)=0$. Therefore, we are considering a time-varying network topology or a temporal network.^{19–21} Considering a temporal network is important in the field of epidemiology because considering temporal information makes a contact analysis of disease spread more accurate.^{26} The non-Markovian structure for edge occurrences may facilitate or depress diffusion;^{27} in particular, interaction patterns in a hospital facilitate the spreading of an epidemic.^{28}

Thus, here, we construct a model that calculates the expected number of infected persons at each day $t$. Our model may be regarded as a discretized version of the Susceptible–Exposed–Infected–Susceptible (SEIS) model,^{12} taking into account an infection history.

In our model, once a person is infected, then the person is called exposed. After $K$ days, the exposed person becomes infectious, being able to infect others. This infectious state will last for $Q\u2212K$ days at maximum and shift to a susceptible state again. Thus, the person becomes susceptible again. We also call people infected if they are either exposed or infectious. For each person, state $si(t)$ at time $t$ corresponds to a $Q$-dimensional vector, where $Q$ is the maximum number of days a person will stay in an infectious state after an infection. According to the statistics of COVID-19, $Q$ seems to be around two weeks.^{6,11} The first component of $si(t)$ corresponds to the first day of the infection, which is an exposed state. The second component of $si(t)$ corresponds to the second day of the infection. Thus, the $q$th component of $si(t)$ corresponds to the $q$th day of the infection for $q=1,2,\u2026,Q$. Let $\sigma k$ be the rate at which a person will keep the exposed or infectious state for $k$ days after the infection. Here, we assume that $0<\sigma k<1$. Then, for person $i$, we can have the following dynamics for a day using a $Q\xd7Q$ matrix $Pi$:

For simplicity, we here assume that $\sigma k$ is constant and $\sigma k=\sigma $. This matrix is similar to the model by Leslie;^{29} thus, we can consider an infection history naturally. The matrix by Leslie was considered in an epidemic problem previously,^{30} especially in an epidemic problem for deers.^{31}

We also use the following $Q\xd7Q$ matrix $\gamma i,j(t)$ for describing the process of infection:

where we have $\gamma i,j(t)(1,k)=0$ for $k=1,2,\u2026,K$, which corresponds to an exposed state, and $\gamma i,j(t)(1,k)=\gamma i,j,k$ for $k=K+1,K+2,\u2026,Q$, which corresponds to an infectious state in the context of the SEIS model. It has been made clear that there is a median interval of 4.6 days between successive infections.^{32} As an infection matrix from person $j$ to person $i$, we have

where $\gamma i,j,k(t)$ is the infection rate specific to susceptible person $i$ at day $t$ while the infectious person $j$ is at $k$ days after the infection. These definitions mean that an exposed person becomes infectious after $K$ days and then the infectious person may infect another susceptible person between $(K+1)$ and $Q$ days after an infection. If exposed, Eq. (2) will contribute to the first component of $si(t+1)$. In addition, this $\gamma i,j,k(t)$ may depend on how the potentially being exposed person $i$ as well as the potentially being infecting person $j$ take their countermeasures for the infection and how long the person $i$ spends after their infection. The countermeasures include, but not limited to, wearing masks, gloves, and/or face shields as well as washing hands and other parts of their bodies.

Now, we are ready to define the whole system. Let $s(t)$ be a vector combining all the states for everyone,

Here, $\u2032$ means the transpose of a matrix. This $s(t)$ is an $NQ$-dimensional vector. Let $T(t)$ be the infection dynamics at day $t$,

This $T(t)$ is an $NQ\xd7NQ$ matrix. Then, our day-resolution equation is as follows:

Here, $s(i\u22121)Q+q(t)$ is the probability that at day $t$, person $i$ is at $q$ days after an infection. From day $t0$, the model can be simulated as

Thus, the total expected number of infected people at day $t$ is written as

where the function $\Gamma $ is defined as

This model is naturally taking into account an infection history of an epidemic, namely, how many days have passed after the infection. In the current setting, we are allowing the recovered people to be exposed more than once. This model may count each infected person more than once if the infectious person has contacts with the other infected persons. However, the total number of the infected persons shown in Eq. (8) does not exceed the number $N$ of persons in this network.

It is easy to derive the upper bound of the dominant eigenvalue for $T(t)$ as

For the derivation of the above upper bound, see the Appendix. From the viewpoint of suppressing the epidemic disease, the upper bound of Eq. (10) should be made as small as possible. This upper bound may be regarded as the upper bound of the basic reproduction number.^{12} In particular, if the upper bound is less than 1, we have

Since we assumed $0<\sigma <1$, this inequality can be simplified as

If Eq. (12) is satisfied, it implies that the basic reproduction number^{12} is less than 1; thus, we can suppress the spreading of the epidemic disease.

If we assume $\gamma i,j,k(t)=\gamma k(t)$, then we have

meaning that

This inequality means that hub persons should have fewer contacts than the theoretical bound $nc$ per day. An important point of Eq. (14) is that the left hand side only depends on contacts by people, while the right hand side is mainly related to the properties of a disease. This $nc$ is the maximum capacity of the number of persons we are allowed to see each day for each person for any initial conditions of an epidemic.

Note that $\u2211k=K+1Q\gamma k(t)=G(t)$, where $G(t)$ is the sum of $\gamma k(t)$ at day $t$. This $G(t)$ can be interpreted as *the total number of infected persons per contact.* Thus, if we define $\gamma k(t)=G(t)fk(t)$, where $fk(t)\u22650$ and $\u2211k=K+1Qfk(t)=1$, then we can change Eq. (14) to

As the units for the left hand side of Eq. (15) become contacts/patients, the units for $G(t)$ are patients/contacts. Thus, if we let $R0$ be the reproduction number (patients/patients) and $M$ be the number of contacts per patient (contacts/patients), then on average, we may have

The excessive number of contacts per person for each day is defined as

We call this quantity the excessive contacts.

Moreover, if $\gamma i,j,k(t)=\gamma (t)$ is constant across all the persons for each time after being infectious, then $nc$ of Eq. (15) becomes

Furthermore, by assuming $\gamma i,j,k(t)=\gamma (t)$, Eq. (14) can be changed as

For example, if $\gamma i,j,k(t)=\gamma (t)$ and the hub persons see at most three persons each day with $Q=14$, $K=0$, and $\sigma =0.9$, the infection rate $\gamma (t)$ should be less than 3.3% for suppressing the epidemic spread. Thus, we can control two things for reducing the spreading of the epidemic disease: we can reduce the number of people to see for each day for each person; we can reduce the infection rate $\gamma (t)$ to allow the larger capacity of persons to see.

If the underlying contact network is scale-free,^{33} $maxj\u2211i\u2260jai,j(t)$ becomes large, and it is likely that we cannot stop spreading of the epidemic disease. If the underlying contact network is a typical small-world network in the sense of Watts and Strogatz^{34} and we have $\u2211i\u2260jai,j(t)$ small as well as almost constant over $i$, we can possibly suppress the spreading of the epidemic disease,^{35} while keeping the minimal interactions to sustain our daily lives.

Therefore, to design countermeasures for stop spreading an epidemic disease as well as sustaining our economic activities, we should introduce constraints on the maximal allowed numbers of people we can see for each day for each person by looking for an optimal contact network for day $t$ that maximally stimulates our economies.

First, we verified the above analytical results with numerical simulations. In our numerical simulations, we first prepared two kinds of *static* networks: for each of the scale-free networks, we followed Barabási and Albert^{33} to generate an underlying network randomly. Thus, the generated networks are intended to be neither assortative nor disassortative.^{4} For small-world networks, we first prepared the one-dimensional regular ring, where neighbors within distance 2 are connected. Then, we replaced a local connection with a long distance connection randomly with the probability $0.05$. Furthermore, we set $N=500$, $Q=14$, and $\gamma i,j,k(t)=0.04$. We assume that the spread of an epidemic starts from a person for each simulation.

When we varied $K=0,1,\u2026,13$ and $\sigma =0.85,0.9$, and $0.95$, the infected persons grew as shown in Figs. 1 and 2. These figures show that the importance for the exposed state but not the infecting state as well as the rate at which the infected state is kept. In particular, if $K$ is larger, the infection is more suppressed, while the infection is more likely to spread if $\sigma $ is larger.

These results can also be summarized as in Fig. 3. Whether an outbreak occurs or not depends on each simulation (Fig. 3). However, the largest eigenvalue of $T(t)$ is correlated well with the number of the infected persons with the probability higher than $0.2$ at day $100$ [Fig. 3(a); correlation coefficient: $0.2925$; Table I]. This correlation is known in the existing literature, for example.^{36} In addition, we found that the number of infected people correlates positively with the maximum number of contacts per person for each day [Fig. 3(b); correlation coefficient: $0.2206$] as well as negatively with the theoretical threshold $nc$ [Fig. 3(c); correlation coefficient: $\u22120.2079$]. As a total, whether or not the excessive contacts [Eq. (17)] are greater becomes a good indicator for forecasting whether or not the number of infected persons at day 100 is 1 or greater [Fig. 3(d); Table II]. Moreover, please observe that Eq. (14) holds numerically as well: If the excessive contacts were negative, then outbreaks did not occur. However, if the excessive contacts become positive, outbreaks may have happened especially certainly for the small-world network examples. There are many cases (5655 cases) in the simulations shown in Table II, where the number of infected people at day 100 is less than 1 even if the excessive contacts are positive. To evaluate such cases more closely, we examined all the 100 cases each of $(\sigma ,K)=(0.9,2)$ and $(0.95,4)$ for the scale-free networks and the 100 cases each of $(\sigma ,K)=(0.9,1)$ and $(0.95,2)$ for the small-world networks using 100 different initial conditions for each network topology. We found that in these networks, the largest eigenvalue of $T(t)$ is greater than $1$ and the excessive contacts are positive. In addition, for each of these networks, there are some initial conditions from which a pandemic starts. Therefore, under the current COVID-19 situation where we do not have much direct treatment options, there is a certain value for a safer guideline, like Eq. (14) we introduced here, that each susceptible person can follow easily for suppressing the epidemic spread for sure.

The maximum for the absolute value of the eigenvalues . | The number of infected people is less than 1 . | The number of infected people is 1 or greater . | Total . |
---|---|---|---|

1 or smaller | 100% (6060) | 0% (0) | 100% (6060) |

Greater than 1 | 36% (849) | 64% (1491) | 100% (2340) |

The maximum for the absolute value of the eigenvalues . | The number of infected people is less than 1 . | The number of infected people is 1 or greater . | Total . |
---|---|---|---|

1 or smaller | 100% (6060) | 0% (0) | 100% (6060) |

Greater than 1 | 36% (849) | 64% (1491) | 100% (2340) |

The excessive contacts . | The number of infected people is less than 1 . | The number of infected people is 1 or greater . | Total . |
---|---|---|---|

Zero or negative | 100% (1254) | 0% (0) | 100% (1254) |

Positive | 79% (5655) | 21% (1491) | 100% (7146) |

The excessive contacts . | The number of infected people is less than 1 . | The number of infected people is 1 or greater . | Total . |
---|---|---|---|

Zero or negative | 100% (1254) | 0% (0) | 100% (1254) |

Positive | 79% (5655) | 21% (1491) | 100% (7146) |

Tables I and II summarize the results of the largest eigenvalue of $T(t)$ and those of excessive contacts, respectively. Both criteria are fine because the outbreaks do not occur if these criteria are less than the corresponding thresholds. On the other hand, comparison between Tables I and II imply that the largest eigenvalue of $T(t)$ provides more reliable information for the outbreaks than the excessive contacts of Eq. (17). However, we need to know all the global information to derive the largest eigenvalue of $T(t)$, while local contact information is sufficient for each individual person to examine the excessive contacts of Eq. (17). Thus, an advantage for using the excessive contacts is that each individual person can easily follow the criterion and decide how they should behave in each day.

Moreover, we also conducted simulations on time-varying small-world networks or, namely, temporal networks (Fig. 4). In these simulations, for each day, we generated a small-world network, which is perturbed from the one-dimensional ring as described above. Here, we set $L=500$, $K=4$, and $\gamma i,j,k(t)=0.04$. Then, we ran the model for the time forward by applying different contact network structures every day. Therefore, this simulation can be regarded as a simulation within which we usually communicate with neighbors, but we sometimes travel for a long distance. In addition, we also assumed that there is a big gathering where half of the population gets together once in 7, 14, 21, and 28 days. The results are shown in Fig. 4. Even in this case, outbreaks of the epidemic were suppressed if the large gathering is once in 4 weeks. This property can be estimated from the following argument based on the means of the maximum eigenvalues for the underlying dynamics. The mean for the maximum eigenvalue of the usual day was about $0.96$, while the mean for the day of the large gathering was about $1.69$. Thus, the maximum eigenvalue for a $\tau $ day cycle of a big event can be estimated as $0.96\tau \u22121\xd71.69$. By solving $0.96\tau \u22121\xd71.69=1$, we found $\tau \u224814.2$. Since the actual maximum eigenvalues fluctuate, this approximation roughly agrees with the results obtained in Fig. 4.

Furthermore, we conducted simulations using a more realistic $\gamma k(t)$ for the COVID-19,^{11} which follows the Weibull distribution^{37} of a shape parameter of $2.826$ and a scale parameter of $5.665$. Then, the shape of infecting probability depending on the number of days after an infection looks as Fig. 5.

We varied $G(t)$ between $0$ and $2$ randomly and ran simulations. We used either scale-free networks or small-world networks. We call the case of scale-free networks as the first type. For each step, we generated a scale-free network randomly as we did above. When we used small-world networks, we prepared two types of simulations, which are the second and third types. In the second type, we generated a network structure for each time starting from a ring, where neighbors within distance 2 are connected. In the third type, we generated a network structure similarly in the second type except that we started from a ring, where only neighbors within distance 1 are connected. In each type, we ran 300 simulations by setting the number of nodes to $500$.

Then, we obtained the results shown in Fig. 6. We found that there are clear epidemic thresholds of $G(t)\u22480.5,0.5,1.0$ for the first, second, and third types of simulations, respectively. These epidemic thresholds of $G(t)$ seem to be inversely related to the mean numbers of newly infecting local directions, which are $2$, $2$, and $1$, respectively. The mean number of newly infecting local directions in the scale-free networks was obtained from the mean number of connections per node, while the corresponding mean numbers for the small work networks were obtained from the initial ring configurations as shown in Fig. 7. In addition, we found in Fig. 8 that the mean number of newly infecting local directions almost agrees with theoretical bounds $nc$ to decide whether an epidemic occurred or not. In particular, the epidemic thresholds for $nc$ for the first, second, and third types of simulations were about $2$, $2$, and $1$ numerically, respectively. Therefore, these simulations show that our theoretical bound $nc$ is a good means for distinguishing whether an epidemic is likely to occur or not. In particular, if the theoretical bound $nc$ is larger than the mean number of newly infecting local directions, then an epidemic is likely to be prevented.

If we follow the above estimation as well as Eq. (16), we may conduct the following calculation. In the early 2020 in China, the reproduction number $R0$ was identified as $2.0$.^{11} The number $M$ of contacts per person in China in 2015 was estimated as $14.37801$ using the contact information of Ref. 38 and the demographic data at the United Nations’ site.^{39} If we assume that the underlying network is scale-free and drawn randomly for each time, $G(t)$ could be $G(t)\u2248R0M=2.0/14.37801\u22480.139$. This value implies that the hub person may be able to meet at most $nc=1/G(t)=1/0.139\u22487.189$ persons per day for suppressing an epidemic. If the underlying network is scale-free, the upper bound should be larger than $7.189$. Although this estimated upper bound should be confirmed with actual data in the future, the number could convey how much we should restrict our daily activities.

We also applied the above statistics obtained from China in the early 2020 and examined whether or not COVID-19 will spread in a hospital based on the real data of human interactions in Ref. 40. We found that the maximum numbers of contacts per person per day in this hospital during the five days were 26, 40, 36, 34, and 34, respectively. Since these numbers are larger than the number of $7.189$ shown above, there is a certain probability that COVID-19 would spread in the hospital if we do not take any countermeasures. Thus, if we have a contact network for a day among persons in a group, we can evaluate whether or not COVID-19 will be likely to spread in the group.

The strength of our model is that we are allowing the uncertainties for the underlying epidemic dynamics including its stochasticity. We only assume on how likely we get exposed when contacted with an infected person, $\gamma (t)$, and how we recover from an epidemic disease, namely, $\sigma $, $K$, and $Q$. Our model can be simulated very realistically by using an actual contact network for each day, considering each individual person. However, we can also work on it theoretically because it is simple enough.

Our criterion of Eq. (14) is a sufficient condition for suppressing epidemic spread for any initial conditions. The meaning is clear because the left hand side depends on contacts by people only, while the right hand side is mainly related to the properties of an epidemic disease. In practice, some length of time is necessary for the epidemic to reach the hub persons before its spreading. In the above simulations of static scale-free networks of $(\sigma ,K)=(0.9,1)$, we have checked that a pandemic did not happen at day 100 if the hub persons are more than five links away from the source person of the infection. We are not sure how much we could generalize this result for a larger network. However, this observation implies how much the hub persons should be protected away from the infection. When we plan countermeasures practically, we do not have to consider all the persons in the whole area but a certain isolated cohort of persons for evaluating Eq. (14). It would be also an interesting open question to discuss the relationship between the microscopic quantity of Eq. (14) and a more macroscopic quantity of the epidemic threshold previously discussed by Speidel *et al.*^{41} on temporal networks.

Our work is not the first paper analytically relating a temporal network with an epidemic (see a recent review^{42} on the topic bridging between epidemics and temporal networks). For example, Huang and Yu^{36} considered a model that can be regarded as a simpler version of the current model; Surano *et al.*^{43} tried to infer strong ties from spreading data; and Bodych *et al.*^{44} formulated a mathematical model, where $k$ contacts are necessary to be infected. The epidemic spread may be facilitated if the older connections are more likely to appear,^{45} while the epidemic spread may be made slower if successive inter-interaction intervals are positively correlated.^{46} However, as far as we have searched, there has been no paper yet that considered an individual-based temporal network model that contains the history of infection and yields a criterion of Eq. (14) or similar, which can be interpreted by each susceptible individual person. This criterion interpretable by each *susceptible* individual person is the point we wanted to highlight in this paper.

In the future work, we will discuss an economic policy based on this model. At this stage, the above simple individual-based criterion becomes the key ingredient for the constraint.

In sum, we have introduced a topological model of epidemic disease spreading. The model is a time-varying linear matrix equation, allowing daily changes for the underlying contact network. By obtaining the upper bound for the dominant eigenvalue of the underlying dynamics, we discuss the implications for the obtained results. We found that hub persons are key for suppressing the spread of an epidemic disease. If we assume that the infection rate is common among all the persons, there is the maximal allowed number of persons whom we can see on each day for each person. Our simple message for the current COVID-19 pandemic is “keep contacts as few as possible.”

## SUPPLEMENTARY MATERIAL

See the supplementary material for the MATLAB codes that are used for generating simulations of Figs. 1–4, 6, and 8 in this study.

We deeply appreciate Ms. Satoko Umino, the conversations with whom have brought us benefits during the manuscript preparation. Our research was supported by the JSPS KAKENHI (Grant No. JP18K11461).

## DATA AVAILABILITY

The data that support the findings of this study are available within the article and its supplementary material.

### APPENDIX: DERIVATION FOR THE UPPER BOUND FOR THE LARGEST EIGENVALUE

Here, we show how to obtain our upper bound of Eq. (11).

Suppose that a non-negative square matrix $B=(bi,j)$ of size $M\xd7M$ is given. Then, the largest eigenvalue for $B$ is upper bounded by $min[maxj\u2211i=1Mbi,j,maxi\u2211j=1Mbi,j]$.

We slightly modify the lines of a similar statement in Ref. 47.

Let $\lambda $ be the largest eigenvalue for $B$ and $x=(xj)$ and $y=(yi)$, the corresponding right and left eigenvectors, respectively. We can generally assume that $x\u22600$ and $y\u22600$. Observe that $x$ and $y$ share the same largest eigenvalue $\lambda $ because their characteristic equations^{48} coincide [i.e., $det(B\u2212zI)=det(B\u2032\u2212zI)=0$].

By using Lemma 1, it is straightforward to obtain Eq. (11).

## REFERENCES

*Proceedings of IEEE 24th Annual Joint Conference of the IEEE Computer and Communications Societies*(IEEE, 2005).

*Gyoretsu No Kaisekigaku (Calculus of Matrices)*(Tokyo Tosho, Tokyo, 2017) (in Japanese).