Modeling Emerging Interpersonal Synchrony and its Related Adaptive Short-Term Affiliation and Long-Term Bonding: A Second-Order Multi-Adaptive Neural Agent Model

1. Introduction

Whenever people interact, their behavior tends to become mutually coordinated in time, or synchronized. Interpersonal synchrony has been found to enhance relationship functioning, for example, by inducing greater levels of closeness, concentration, coordination, cooperation, affiliation, alliance, connection, or bonding.^{1,2,3,4,5,6,7,8,9,10} In literature^11,12 it was suggested that to model the complex, cyclical types of dynamics that occur, a dynamical systems modeling approach is needed. Notably, the benefits of interpersonal synchrony include patterns of mutual adaptation both in the short term and in the long term. For instance, in the context of psychotherapy, a patient and therapist who synchronize their movements⁸ may experience a stronger sense of sharing the present moment during a therapeutic session.¹³ Over multiple sessions, this increased social presence may strengthen the therapeutic bond, which allows the patient and therapist to work together more effectively.¹⁴

This paper has as main goal — to address computational and mathematical analyses of the complex adaptive dynamics of such forms of short-term and long-term adaptivity of interaction behavior related to interpersonal synchronization and to verify the hypothesis that the underlying mechanisms put forward in the literature indeed generate these social phenomena by an emerging and adaptive interactive process. More specifically, these analyses cover three different but closely related levels:

•	Analysis of mechanisms from the literature in psychology and neuroscience that are suggested to play a role in these complex multi-adaptive dynamics. These include mechanisms for different forms of (synaptic and nonsynaptic) plasticity and control over plasticity by metaplasticity (also called second-order plasticity or second-order adaptivity). Here, the underlying hypothesis in the literature is that these mechanisms put forward are sufficient and able to generate the emerging and adaptive patterns of synchronization and adaptation of the interaction behavior. This hypothesis is tested in silico in this paper by computational simulation based on these mechanisms.
•	Analysis of mathematical formalization of these mechanisms behind the considered complex adaptive dynamics by an agent-based second-order multi-adaptive dynamical system. This covers analysis of conducted simulation experiments based on such formalization and includes analysis of stationary points and equilibria for the occurring dynamics and adaptivity.
•	Analysis of how the specific representation of such an adaptive dynamical system by self-modeling networks used here is positioned in the wider landscape of adaptive dynamical systems. It is analyzed how the self-modeling network format can be used to provide a canonical representation for any smooth adaptive dynamical system, which also covers most neural system models.

So, more specifically, the neural agent model that is a central focus here is an adaptive dynamical systems model based on a number of mechanisms in the literatures on cognitive, behavioral, and affective neuroscience. A neural basis for short-term behavioral adaptivity can be found in the recent work on the (nonsynaptic, intrinsic) adaptive excitability of (neural) states.^15,16,17,18 By contrast, a neural basis for long-term adaptivity can be found in the classical notion of synaptic plasticity.^19,20,21,22 Together, these two fundamentally different forms of adaptation yield a model of a multi-adaptive neural agent. The two forms of adaptation also interact with each other.

The extent of adaptation that an agent requires may vary from situation to situation. The capacity to adjust plasticity to the demands of the situation relates to metaplasticity.^23,24 This model of a neural agent models metaplasticity as a second-order form of plasticity that controls plasticity in a context-sensitive manner. The resulting model yields a second-order multi-adaptive neural agent, which is human-like in the sense that it incorporates an interplay of three major mechanisms for adaptivity that according to the neuroscientific literature characterizes human agents.

Note that it is not claimed that the model is human-like in the sense that it would cover the applied neural mechanisms at a physiological level addressing that level was not within the scope of this research. Instead, these neural mechanisms were considered and modeled in a more abstract manner at a functional level. Investigating physiological scalability would be another, next enterprise to be addressed.

This model of a neural agent further includes intrapersonal synchrony and interpersonal synchrony and their links to short-term and long-term behavioral adaptivity. To model the pathway from synchrony patterns to this behavioral adaptivity, we included both built-in intrapersonal synchrony and interpersonal synchrony detectors. Here, intrapersonal synchrony means that within an agent, actions for the different modalities occur in a coordinated manner. Interpersonal synchrony means that for each modality, the actions of the two agents occur in a coordinated manner. The addressed modalities are movement, affect, and verbal modalities. We included these three modalities because they have each been shown to be influential in interpersonal behavior.¹⁴

We evaluated the neural agent model in a series of simulation experiments for two agents with a setup in which a number of stochastic circumstances were covered in different (time) episodes. The simulations included not only episodes with a stochastic common stimulus for the two agents, but also episodes with different stochastic stimuli for the agents. Moreover, to analyze the role of communication, stochastic circumstances were also included for episodes when communication was enabled by the environment and episodes when communication was not enabled.

Next, as part of further analysis of the self-modeling network modeling approach it is shown how any (smooth) adaptive dynamical system can be modeled in a canonical way as a self-modeling network model. In this way, any adaptive dynamical system has its canonical representation as a self-modeling network model and can be analyzed based on this canonical representation. On the one hand, this shows that the chosen modeling approach does not introduce biases or limitations if adaptive dynamical systems are modeled using it. In particular, it shows also that the approach generalizes most common neural system models. On the other hand, this was a basis to show how stationary point analysis and equilibrium analysis for adaptive dynamical systems can be performed by using the self-modeling network representation for the specific adaptive dynamical system model introduced here.

2. Main Assumptions and Background Knowledge

In this section, we present the main assumptions behind the introduced adaptive neural agent model and relate them to the relevant neuroscience literatures. This grounding in neuroscience is based on pathways for a circular interplay of synchrony with both nonsynaptic plasticity¹⁶ and synaptic plasticity,^19,20,21,22 thereby covering both short-term time scales and long-term time scales and their interaction. More specifically, the following underlying assumptions are made for the pathways involved; for a conceptual overview, see Fig. 1. Note that the main example used in the presentation concerns two agents A and B and their interaction.

2.1. Interpersonal synchrony leads to adaptation of interaction behavior

Interpersonal synchrony is often followed by a behavioral change or adaptation of mutual behavior.^{1,2,3,4,5,6,7,8,9,10} This adaptive shift in mutual behavioral coordination has been observed, for instance, in psychotherapy sessions. Research has shown that therapists were rated more favorably and as more empathic when, beforehand, they were instructed to make their movements more synchronized with the client.^25,26,27 Similarly, Ramseyer and Tschacher⁸ found that initial movement synchrony between client and therapist was predictive of the client’s experience of the quality of the alliance at the end of each session. Furthermore, Koole and Tschacher¹⁴ reviewed converging evidence that movement synchrony has a positive effect on the working alliance between patient and therapist. More generally, synchrony in face-to-face interactions has been found to promote interpersonal affiliation.^10,28

2.2. Behavioral adaptation after interpersonal synchrony occurs both in the form of short-term adaptation and long-term adaptation

Much research on interpersonal synchrony has focused on short-term adaptive changes in interpersonal coordination.^1,6,9,10,29 However, several lines of research have observed effects of interpersonal synchrony on long-term adaptation as well. First, developmental research has observed that movement synchrony between infant and caregivers predicts social interaction patterns of the child several years later.²⁸ Second, research on close relationships suggests that early patterns of interpersonal synchrony predict subsequent indicators of relationship functioning, for instance, one study found that spouses’ patterns of cortisol variation converged over a period of years, indicating long-term shifts in interpersonal coordination.³⁰

Third and last, research on psychotherapy processes has found that markers of interpersonal synchrony in early sessions can predict the development of the therapeutic relationship⁸ and therapeutic outcomes.⁷ Long-term adaptation processes remain less well-studied than short-term adaptation processes. Nevertheless, the convergence of evidence is sufficient to conclude that interpersonal synchrony is likely to promote both short-term and long-term adaptation in interpersonal relationships.

2.3. The behavioral adaptation relies on different neural mechanisms: Synaptic plasticity of connections and nonsynaptic plasticity of intrinsic excitability

In the neuroscientific literature, a distinction is made between synaptic and nonsynaptic (intrinsic) adaptation. The classical notion of synaptic plasticity has been used to explain long-term behavioral adaptation.^19,20,21,22 This addresses how the strength of a connection between different states is adapted over time due to simultaneous activation of the connected states. By contrast, the nonsynaptic adaptation of intrinsic excitability of (neural) states has been addressed in more detail more recently.^15,16,17,18 The latter form of adaptation has been related, for example, to homeostatic regulation¹⁷ and also to how deviant dopamine levels during sleep make that dreams can use more associations due to easier excitable neurons.³¹ Moreover, both (synaptic and nonsynaptic) forms of adaptation can easily work together.³²

In the neural agent model these two adaptation mechanisms and their interaction have been used to model behavioral adaptivity: the former for long-term adaptation and the latter for short-term adaptation. Here an interplay of two types of adaptivity occurs. Synchrony does not only lead to short-term adaptation, but short-term adaptation itself also intensifies interaction which can lead to more synchrony which in turn can strengthen the long-term adaptation. Besides, long-term adaptivity also strengthens interaction which leads to more synchrony and consequently stronger short-term adaptivity. In this way, via multiple circular pathways a dynamic interplay occurs between synchrony, short-term adaptivity and long-term adaptivity.

Plasticity is not a constant feature, as it often is highly context-dependent according to what is called metaplasticity.^23,24 For example, “adaptation accelerates with increasing stimulus exposure”.²⁴ Toenable such context-sensitive control of plasticity, second-order adaptation (i.e. adaptation of the adaptation) has been included in the neural agent model, which makes the model more realistic.

2.4. The pathways from synchrony to behavioral adaptation involve synchrony detection states

If synchrony occurs for an agent and due to this the agent adapts the interaction behavior, this suggests that agents possess a facility to notice or experience synchrony patterns for the different modalities. Indeed, the assumption is made that agents do in some way (perhaps unconsciously) detect synchrony and from there may trigger behavioral adaptation for their interaction behavior. In the pathway from synchrony patterns to changed interaction behavior patterns, such synchrony detection states can be considered as specific mediating mental states. Such a state p in general has been called a mediating state for the effect of a past pattern a on a future pattern b entailed by pattern a (Refs. 33 and 34); similarly, such a (brain) state is referred to as describing “informational criteria” for future activation.^35,36 In line with previous research,³⁷ it is assumed that not only the detected interpersonal synchrony but also the detected intrapersonal synchrony relating to a conscious emotion has a causal effect on the behavioral adaptivity.

3. Self-Modeling Network Modeling

The presented neural agent model is an adaptive dynamical system model designed and specified based on network-oriented modeling. The network-oriented modeling approach used here is basically a causal network modeling approach where nodes model states that have activation values that change over time and connections between these states model causal relations that have their effects in a temporal, dynamic manner on the state activations. Thus, dynamical systems are modeled. Moreover, by enabling these causal relations and characteristics for their effects on state activations to change as well, also adaptive dynamical systems are covered. This is also done in a network-oriented manner, using a so-called self-modeling network architecture. In Sec. 6, it will be shown that any smooth adaptive dynamical system can be modeled in this way. In this section, this modeling approach is briefly introduced.

Following the network-oriented modeling approach^38,39,40 used here, a temporal–causal network model is characterized by (here X and Y denote nodes of the network, also called states, and $X(t)$ and $Y(t)$ denote their activation values at time t):

•	Connectivity characteristics Connections from a state X to a state Y and their weights ${\omega }_{X,Y}$.
•	Aggregation characteristics For any state Y, some combination function $c_{\pi Y,Y}$ ($V_1,\dots ,V_k$) with vector of parameter values ${\pi }_Y=({\pi }_{1,Y},\dots ,{\pi }_{m,Y})$ defines the aggregation that is applied to the impacts $V_i={\omega }_{X_i,Y}{X}_i(t)$ on Y from its incoming connections from states $X_i$.
•	Timing characteristics Each state Y has a speed factor ${\eta }_Y$ defining how fast it changes for given causal impact.

Note that for the sake of notational simplicity, often $c_{{\pi }_Y,Y}$ will be denoted by $c_Y$, omitting the subscript ${\pi }_Y$; this will not mean that there are no parameters, they are just left implicit.

These network characteristics ${\omega }_{X,Y}$, $c_{{\pi }_Y,Y}$, ${\pi }_Y$, and ${\eta }_Y$ for a given network model serve as a (formal) design specification of this network model. The following canonical difference and differential equations for temporal–causal network models used for simulation and analysis of such network models, incorporate these network characteristics ${\omega }_{X,Y}$, $c_Y$, ${\pi }_{i,Y}$, and ${\eta }_Y$ in a standard numerical format :

for any state Y and where $X_1$ to $X_k$ are the states from which Y gets its incoming connections. Note that (1) has a format similar to that of recurrent neural networks. Within the dedicated software environment implemented in MATLAB, a large number of currently around 60 useful basic combination functions are included in a combination function library.³⁹

The above concepts enable us to design network models and their dynamics in a declarative manner, based on mathematically defined functions and relations and specified in a standard table format covering all network characteristics (called role matrices, see Appendix  B). The examples of combination functions that are applied in the model introduced here can be found in Table 1. Here, for the third and fourth function, rand(1, 1) draws a random number from [0, 1] in a uniform manner and a is a persistence factor (with value 0.5 used in the simulations).

Realistic network models are usually adaptive: often not only their states but also some of their network characteristics change over time. By using a self-modeling network (also called a reified network), a similar network-oriented conceptualization can also be applied to adaptive networks to obtain a declarative description using mathematically defined functions and relations for them as well.^39,40 This works through the addition of new states to the network (called self-model states) which represent (adaptive) network characteristics. In the graphical three-dimensional (3D)-format as shown in Sec. 4, such additional states are depicted at a next level (called self-model level or reification level), where the original network is at the base level.

As an example, the weight ${\omega }_{X,Y}$ of a connection from state X to state Y can be represented (at a next self-model level) by a (connectivity) self-model state named $W_{X,Y}$, which can be used to model synaptic plasticity.^19,20,21,22 Similarly, all other network characteristics from ${\omega }_{X,Y}$, $c_Y$ (…) and ${\eta }_Y$ can be made adaptive by including self-model states for them. For example, for adaptive excitability^15,16,17 the threshold ${\tau }_Y$ (for a logistic combination function) for a state Y can be represented by a (aggregation) self-model state named $T_Y$ and an adaptive speed factor ${\eta }_Y$ can be represented by a (timing) self-model state named $H_Y$. Dynamics for the activation values for these self-model states are modeled by adding their own network characteristics (thus integrating them in the network structure) and applying Eq. (1) for them.

If for all network characteristics $\omega$, $\pi$, $\eta$ for all base level states, respective self-model states W, P, H are introduced representing these network characteristics, then the canonical difference and differential equations for the base level states of the self-modeling network model are

where $P_Y(t)=(P_{1,Y}(t),\dots ,P_{m,Y}(t))$.

This canonical difference equation is incorporated in the dedicated software environment. By instantiating this general difference equation (2) by proper values for the network characteristics for all base states Y and similarly instantiating equation (1) for all self-model states, the software environment runs a system of n difference equations where n is the number of (base and self-model) states in the network.

When Eqs. (2) are compared to Eqs. (1), it can be noticed that at each point in time t, for the value of each network characteristic the activation value of its corresponding self-model state is used: for ${\eta }_Y$ the value $H_Y(t)$ is used, for ${\omega }_{X_i,Y}$ the value $W_{X_i,Y}(t)$, etc.In this way, each of these self-model states is assigned the functional role of the specific network characteristic it represents. In particular, when the activation values of these W-states, P-states, and H-states change, the corresponding network characteristics change accordingly. This makes these network characteristics adaptive.

More mathematical background of this self-modeling network architecture construction and Eqs. (2) has been described elsewhere.³⁹

Usually not for all network characteristics self-model states are introduced but only for part of them. For example, in case only self-model states P for the combination function parameters $\pi$ are introduced, then the canonical difference and differential equation are

where $P_Y(t)=(P_{1,Y}(t),\dots ,P_{m,Y}(t))$. This specific case will come back in the mathematical analysis addressed in Sec. 6.2 (to establish Theorem 2 there).

Note that difference and differential equations (2) are not exactly in the standard format of a temporal–causal network, as $H_Y$ is not a constant speed factor and also the P-and W-values are not constant. However, it can be rewritten into the temporal–causal network format when the following combination function $c_Y^{*}(..)$ is defined:

where $P=(P_1,\dots ,P_m)$.

Based on this combination function, consider the following difference equation:

This is indeed in temporal–causal network format (1) (with speed factor 1). Now note that using (4), Eq. (5) can be rewritten as follows :

where $P_Y(t)=(P_{1,Y}(t),\dots ,P_{m,Y}(t))$.

Equation (6) shows exactly difference equation (2) above; this confirms that the chosen combination function $c_Y^{*}(..)$ in (4) to show that the self-modeling network has a temporal–causal network format (1) works.

As the outcome of a process of network reification is also a temporal–causal network model itself, as has been shown above, this self-modeling network construction can easily be applied iteratively to obtain multiple orders of self-models at multiple (first-order, second-order, etc.) self-model levels. For example, a second-order self-model may include a second-order (timing) self-model state $H_{W_{X,Y}}$ representing the speed factor $\eta W_{X,Y}$ for the dynamics of first-order self-model state $W_{X,Y}$ which in turn represents the adaptation of connection weight ${\omega }_{X,Y}$. Similarly, a second-order self-model may include a second-order (timing) self-model state $H_{T_Y}$ representing the speed factor ${\eta }_{T_Y}$ for the dynamics of first-order self-model state $T_Y$ which in turn represents the adaptation of excitability threshold ${\tau }_Y$ for Y.

In this paper, this multi-level self-modeling network modeling perspective will be applied to obtain a second-order adaptive network architecture addressing controlled behavioral adaptation induced by detected synchrony. In this self-modeling network architecture, the first-order self-model models the adaptation of the base level network that models behavior, and the second-order self-model level the control over this adaptation. As an example, the control level can be used to make the adaptation speed context-sensitive as addressed by metaplasticity literature.^23,24 For instance, the metaplasticity principle “adaptation accelerates with increasing stimulus exposure”²⁴ formulated by Robinson et al. can easily be modeled by using second-order self-model states; this actually has been done for the introduced model, as will be discussed in Sec. 4.

4. The Adaptive Neural Agent Model

In this section, our adaptive neural agent model is explained in some detail. The controlled adaptive agent design uses a self-modeling network architecture of three levels as discussed in Sec. 3: a base level, a first-order self-model level, and a second-order self-model level. Here the (middle) first-order self-model level models how connections and excitability thresholds of the base level are adapted over time, and the (upper) second-order self-model level models the context-sensitive control over the adaptations. Appendix B provides explanations for all of its states and a full specification of the model.

4.1. Base level

Figure 2 shows a graphic overview of the base level. For each agent, interaction states were modeled: states involved in sensing (indicated by sense) are on the left-hand side of each box, and states involved in execution or expression of actions (move, exp_affect, talk) on the right-hand side. In between these interaction states, within a box are the agent’s internal mental states; outside the boxes are the world states. Note that we assume that each agent also senses its own actions, modeled by the arrows from right to left outside the box.

For each agent, we modeled a few internal mental states such as sensory representation states (rep) and preparation states (prep) for each of the three modalities: movement m, expression of affect b, and verbal action v.

Furthermore, each agent has a conscious emotion state for affective response b (cons_emotion). Each of the mentioned states is depicted in Fig. 2 by a light pink circle shape. For each modality, its representation state has an outgoing (response) connection to the corresponding preparation state and it has an incoming (prediction) connection back from the preparation state to model internal mental simulation.^41,42

Finally, there are the six synchrony detector states (depicted in Fig. 2 by the darker pink diamond shapes) which are introduced here. As in previous research³⁷ we cover three intrapersonal synchrony detection states for the three pairs of the three modalities:

• movement–emotion ($m-b$),

• movement–verbal action ($m-v$),

• semotion–verbal action ($b-v$).

These intrapersonal synchrony detection states have incoming connections from the two execution states for the modalities they address. The conscious emotion state is triggered by incoming connections from the preparation state for affective response b together with the three intrapersonal synchrony detection states.⁴³ In addition, the conscious emotion state has an incoming connection from the verbal action execution state (for noticing the emotion in the verbal utterance) and an outgoing connection to the preparation of the verbal action (for emotion integration in the verbal action preparation).

There are three interpersonal synchrony detection states for the three modalities m, b, and v. Each of them has two incoming connections: from the sensing state (representing the action of the other agent) and the execution state (representing the own action) of the modality addressed.

For a few states and connections, their excitability and connection weights are adaptive depending on detected synchrony: detected synchrony leads to becoming more sensitive to sensing an agent and expressing to that agent (short-term effect) and to connecting stronger to the agent (long-term effect). Here, two different time scales for the adaptations are considered:

•	On the short term, enhancing the excitability of such internal states, so that they become more responsive or sensitive (a form of instantaneous homeostatic regulation).
•	On the long term, making the weights of such connections stronger so that propagation between states is strengthened (a form of a more endurable bonding).

This applies to two types of states and four types of connections in particular, all playing an important role in the interaction behavior of the two agents:

•$$Short-term adaptive excitability for internal states

$$∘$$The representation states for each of the three modalities.

$$∘$$The execution states for each of the three modalities.

•$$Long-term adaptive internal and external connections

$$∘$$The (representing) connections from sensing to representation states for each of the three modalities.

$$∘$$The (executing) connections from preparation to execution states for each of the three modalities.

$$∘$$The (observing) connections from world states to sensing states.

$$∘$$The (effectuating) connections from execution states to world states.

Thus, more synchrony detected will lead to enhanced excitability for these types of states (short-term adaptation) and for these connections to become stronger (long-term adaptation); each type of all these adaptations contributes in its own way (and time scale) to the interaction behavior of the agents. In the short term, more sensitive states for representations will lead to gaining better images of the modalities of the other agent; this will make the sensed signals better available and accessible for the agent. More sensitive states for execution will lead to better expressed own modalities, so that the other agent can sense them better.

Over time and repeated interactions, a stronger (external) observing connection will lead to sensing the other agent better (e.g. turning sensors in the right direction and bending or getting closer to the other), and a stronger representing connection again (but now in a more endurable manner) will make the sensed signals better available and accessible for the agent. Conversely, a stronger executing connection will also contribute (in an endurable manner) to stronger expression and acting toward the other and a stronger effectuating connection to better availability (for the other) of the action effects in the world (e.g. more visible, better hearable by directing and positioning in the right direction, and bending or getting closer to the other). In Sec. 4.2, we discuss in more detail how we modeled these forms of adaptivity and their control using the principle of self-modeling of the network model.

Finally, at the base level some world states are modeled for stimuli s that are sensed by the agents. In the simulations, they have stochastic activation levels. In some episodes one common stimulus is observed by both agents (for example when they physically meet and therefore are in the same environment), but in other episodes the agents receive different stimuli. Furthermore, also the world situation’s suitability for enabling communication between the two agents is modeled by similar stochastic fluctuations. Moreover, two context states are included to model the conditions to maintain excitability thresholds well.

4.2. Modeling adaptation and its control

We modeled adaptation and its control needed in the neural agent model using a “self-modeling network”^39,40; see also Sec. 3. Following what has been described in Sec. 4.1, for a number of states Y adaptive excitability has been modeled via the excitability threshold ${\tau }_Y$ of the logistic function used for these states (see Table 1). Moreover, the strengthening of connections from X to Y has been modeled via adaptive connection weights ${\omega }_{X,Y}$. Following Sec. 3, these adaptations have been modeled in particular through self-modeling for these ${\tau }_Y$ and ${\omega }_{X,Y}$ by adding the following first-and second-order self-model states:

•	First-order self-model T-states TY are used for short-term adaptation of the adaptive base excitability thresholds ${\tau }_Y$ for the internal representation states and execution states Y for the three considered modalities (movement, affective response, and verbal action). For each agent there are six of these T-states, both for the three representation base states and the three execution base states (both for the three modalities).
•	First-order self-model W-states W${}_{X,Y}$ are used for adaptation of the adaptive base connection weights ${\omega }_{X,Y}$ for both internal and external connections for the three considered modalities; internal connections at the base level from sense states to representation states and from preparation states to execution states, and external connections from execution states to world states and from world states to sense states. For each agent there are 12 of these W-states, for the connections from world states to sensor states, from sensor states to representation states, from representation states to execution states, and from execution states to world states (all for the three modalities).
•	Second-order self-model H${}_T$-states are used for control of the T-states for adaptation of the adaptive excitability thresholds ${\tau }_Y$ for the internal representation states and execution states Y. For each agent there is one of these states.
•	Second-order self-model H${}_W$-states are used for control of the W-states for the adaptation of the adaptive base connection weights ${\omega }_{X,Y}$. For each agent there is one of these states.

Figure 3 shows the overall design of the network model; here, the first-order self-model states are in the middle (blue) plane and the second-order self-model states in the upper (purple) plane. The first-order states include T-states representing the excitability thresholds of representation and execution states and W-states representing the weights of the different types of adaptive connections addressed.

Recall from Sec. 3 the canonical difference and differential equation (2) for a self-modeling network. In this equation it is shown that at each time point for the values of the network characteristics the values of these self-model T-states and W-states are used. Based on this equation, by changing the activation values of these T-states and W-states over time t, the corresponding excitability thresholds and connection weights change accordingly which makes them adaptive. Such change of the values of the T-states and W-states occurs due to the influences from the detected synchronies, modeled by the upward (blue) arrows in Fig. 3 from the synchrony detection states in the base plane to the T-states and W-states in the middle plane.

For most states the combination function alogistic is used which has an excitability threshold parameter that can be made adaptive; see the last column in Table 1. The synchrony detections states, however, have a different function called compdiff to measure the extent of synchrony. For further details, see Table 1 and Appendix  B.

There are four second-order self-model states to control the adaptation: two second-order self-model states H${}_{T_A}$ and H${}_{T_B}$ for excitability adaptation control, one for each agent, and two second-order self-model states H${}_{W_A}$ and H${}_{W_B}$ for connection weight adaptation control, also one for each agent. These second-order self-model states are used to represent the adaptation speed (learning rate) for the adaptive excitability threshold T-states and connection weight W-states for the concerning agents A and B. Based on the canonical difference and differential equation (2) for a self-modeling network from Sec. 3, for each time point t the activation values of the second-order self-model states H${}_{T_A}$, H${}_{T_B}$, H${}_{W_A}$, and H${}_{W_B}$ at t are used as the values for these network characteristics of the first-order self-model T-states and W-states for A and B. These second-order self-model states H${}_{T_A}$, H${}_{T_B}$, H${}_{W_A}$, and H${}_{W_B}$ model the second-order adaptation (or metaplasticity) principle “adaptation accelerates with stimulus exposure”.²⁴ To this end they have incoming connections (blue arrows from base plane to upper plane) from the stimulus representation states at the base level.

5. Simulation Results

Appendix B provides a full specification of the model as used in our simulations. As can be seen, in general the values have been chosen in a standard manner. For example, all positive connection weights are 1, except for the long-term adaptation speed self-model states H${}_{WA}$ and H${}_{WB}$ which are 0.01, see Table B.4. The negative connection weights for the T-states are −0.12 so that they together add up to a value −0.72, which compensates the positive weight 1 enough to get effect. The values for the steepness parameter $\sigma$ (see Table B.6) for the combination function alogistic were all set on 5 which also is a kind of default value. For values of the threshold parameter $\tau$, they often have been set at 0.5 but accordingly higher for states with multiple incoming connections.

For all states, when the Euclidean combination function eucl is used, the order parameter n is 1 which makes it linear and the scaling factor $\lambda$ is the sum of the weights of the incoming connections which normalizes it. Note that the time unit is kept abstract. Depending on the application context one might think of minutes, for example for therapeutical sessions.

5.1. Design of the simulation experiments

In this section, we evaluate our neural agent model in an experimental simulation paradigm. Our paradigm was set up in such a way that we could evaluate the behavior of our two agents during four different types of consecutive episodes (see Table 2 and Fig. 4) which are explained as follows. Each of these types of episodes lasts for 30 time units, so that a cycle of four episodes equals 120 time units. Our total simulation run had a duration of 840 time units and the step size ($\mathrm{\Delta }t$) was 0.5, resulting in 1680 computational steps in total for each simulation run. This means that each cycle of four episodes was repeated seven times in each simulation. As it concerns a partly stochastic simulation, we ran 20 repetitions of each simulation with the same episodic paradigm and parameter settings, to get a sense of the robustness of the neural agent model’s behavior. It turned out that general patterns were approximately similar across all independent simulations. Therefore, we selected one simulation to discuss in the upcoming subsections.

Regarding the four different types of episodes in this simulation, they manipulate both whether or not the two agents received the same or a different stochastic stimulus and whether or not they were able to communicate (with some stochastic variations in enabling conditions, due to environmental changes and noise) with each other (Table 2). The specific episodes for the considered example simulation are shown in Fig. 4.

The world states $w{s}_{s,A}$ and $w{s}_{s,B}$ indicate the different stimuli for agent A and B from the world (activated from time 0 to 60 and then repeated every 120 time units; see the dark solid and dashed blue lines for A, respectively B). Similarly, world state $w{s}_s$ indicates the common stimulus (activated from time 60 to 120 and then repeated every 120 time units; see the purple line). These three states have values stochastically fluctuating approximately between 0.7 and 0.9. Furthermore, the self-model states $W_{{\mathrm{\text{exec-ws}}}_{x,A,B}}$ (from A to B) and the states $W_{{\mathrm{\text{exec-ws}}}_{x,B,A}}$ (from B to A) indicate the communication-enabling conditions in the environment. They are activated from time 30 to time 60 thereby fluctuating stochastically roughly between 0.45 and 0.65 and then repeated every 60 time units.

All these stochastic activation patterns indeed follow the pattern shown in Table 2 with repetition every 120 time units until end time 840.

5.2. Behavior of the base states of the neural agent model

For the base states, in the first phase for time 0 to 10 the representations (states rep${}_{s,A}$ and rep${}_{s,B}$) for the stimulus are activated (the curves fluctuating around 0.8) and preparations (states prep${}_{x,A}$) for actions are triggered (curves going to 1); see the upper graph in Fig. 5. This leads, together with the intrapersonal synchrony detection activation (see Figs. 5 and 6), to the conscious emotion around time 10 (red curve going to 1), but this still is only internal processing as no executions of actions take place yet. The action executions (states move${}_{m,A}$, exp_affect${}_{b,A}$, and talk${}_{A,B,v}$) for both agents start to come up after time 10 (e.g. the purple line); this also depends on the short-term adaptations that will be discussed in Sec. 5.4. The curves immediately under these executions concern the sensing of the other agent’s actions (the sense${}_{A,x,B}$ and sense${}_{B,x,A}$ states); in some periods they are slightly fluctuating due to environmental noise on the communication channels. The actual communication level (the ws${}_{x,A,B}$ and ws${}_{x,B,A}$ states) is seen below it from 30 to 60 and from 90 to 120.

For the longer term, the lower graph in Fig. 5 shows that each interval with enabling conditions for communication leads to higher activations of the action executions (the purple line) until values around 0.8 are reached. This is due to a long-term behavioral adaptation that is discussed in Sec. 5.5. Accordingly, the sensing states become higher as well over this longer term, but not as high as the action executions, due to a communication bias incorporated in the model. This overall pattern shows that the enabling conditions for communication have a stronger adaptive effect on the actions than having a common stimulus.

5.3. Behavior of the intrapersonal synchrony and interpersonal synchrony detector states

The curves that the graphs in Figs. 6 and 7 have in common depict the detected intrapersonal synchrony and interpersonal synchrony. Here:

•	The detected intrapersonal synchrony detection is represented by the states intrasyncdet${}_{A,x-y}$ and intrasyncdet${}_{B,x-y}$ shown as the light green and light blue curves going to 1 from time 0 to 15.
•	The interpersonal synchrony detection is represented by the states intersyncdet${}_{A,B,x}$ and intersyncdet${}_{B,A,x}$ shown as the red and blue curves going to 0.4 from time 0 to 30 and further to 0.8 from time 30 to 60.

Here it can be observed that the detection of intrapersonal synchrony takes place already in the first episode from time 0 to 30, meaning no common stimulus or communication is required. In contrast, the detection of interpersonal synchrony strongly depends on the interaction between the two agents. Note also that the former type of detected synchrony reaches a perfect level of 1, due to the coherent internal makeup of the agents, while the latter type does not get higher than around 0.8. At first sight this may look strange, given that the actual executions of actions of both agents are practically the same, as discussed above (see Sec. 5.2 and Fig. 5). However, this is due to the communication bias that was also noted above in Sec. 5.2. This demonstrates the capability of the model that it is able to distinguish a subjective personally detected interpersonal synchrony from an objective form of interpersonal synchrony detection as might be assigned by an external observer but not by the agent itself.

5.4. The interplay between synchrony and short-term adaptation

In Fig. 6, the synchrony detection states are shown together with the states involved in the short-term adaptation: the first-order self-model T-states that represent the adaptive excitability thresholds for representation and execution states and the second-order self-model H${}_T$-states that represent the T-states’ speed factors (adaptive learning rates). Except the synchrony detection states already discussed in Sec. 5.3, the graphs show two light green curves fluctuating around 0.6 for the H${}_T$-states and a blue curve going down to below 0.3 for the T-states.

According to the metaplasticity principle “adaptation accelerates with increasing stimulus exposure”,²⁴ the H${}_T$-states indeed fluctuate with the stimuli. Furthermore, in accordance with this, when one stimulus period is in transition to another stimulus period, it can be seen that there is a short dip in the values of the H${}_T$-states, as stimuli start from 0, so there is a very short period of a lower level, as can also be seen in Fig. 4. Moreover, it is clear that the T-states (e.g. the blue curve) show an opposite pattern compared to the pattern of the interpersonal synchrony detection states. In particular, in the episodes from 30 to 60 and from 90 to 120 (and so on), where the detected interpersonal synchrony is the highest, the T-states for the excitability thresholds are the lowest. This is a short-term adaptation that makes that these agent states related to the communication with the other agent have a higher excitability due to the detected interpersonal synchrony, which will have an intensifying effect on their communication. Not coincidentally, the mentioned periods are also the periods with good enabling conditions for communication (see also Sec. 5.2). It can also be noted that this tendency is a short-term effect and is reversible: the T-states get higher again when the detected interpersonal synchrony gets lower.

5.5. The interplay between synchrony and long-term adaptation

In Fig. 7 the synchrony detection states are shown together with the states involved in the long-term adaptation: the first-order self-model W-states that represent the adaptive weights for the connections to the representation and execution states and the second-order self-model H${}_W$-states that represent the W-states’ speed factors (adaptive learning rates). Here, except the synchrony detection states already discussed in Sec. 5.3, the graphs show the W-states (e.g. a blue curve) slowly and gradually going up to above 0.5 at time 120 and further up to about 0.8 at time 840.

Moreover, at a very low level, the curves for the H${}_W$-states can be seen. They also fluctuate according to the metaplasticity principle “adaptation accelerates with increasing stimulus exposure”,²⁴ but at a very low level around 0.005 (see the lower graph in Fig. 7). Again, following the same principle, when one stimulus period is in transition to another stimulus period, it can be seen that there is a short dip in the values of the H${}_W$-states. This happens because stimuli start from 0, so there is a very short period of a lower level (see Fig. 4). The pattern of the W-states indeed shows a long-term adaptation effect. It highlights that they get a repeated boost in the time intervals 30–60, 90–120, and so on, and show a form of persistency. These boosts occur specifically in these intervals for a reason. These intervals are when there are communication-enabling conditions and as discussed in Sec. 5.4 that induces synchrony and the short-term adaptation via the T-states, which in turn add to synchrony. Therefore, these two effects are at the basis of these boosts for the long-term adaptation. In this way, there is a form of interaction between short-term and long-term adaptation.

6. Modeling and Analysis of Adaptive Dynamical Systems via their Canonical Self-Modeling Network Representation

This section discusses how any smooth adaptive dynamical system can be modeled by a self-modeling network model. It is shown in particular that any adaptive dynamical system has a canonical representation as a self-modeling network defined by network characteristics for connectivity, aggregation, and timing. The network concepts of this canonical representation of an adaptive dynamical system provide useful tools for formal analysis of the dynamics of the adaptive dynamical system addressed. With this idea in mind, equilibrium analysis of self-modeling network models is addressed. Dynamics in network models are described by node states that change over time (for example, for individuals’ opinions, intentions, emotions, beliefs, etc.). Such dynamics depend on network characteristics for the connectivity between nodes, the aggregation of impacts from different nodes on a given node, and the timing of the node activation updates.^39,40

For example, whether within a well-connected group in the end a common opinion, intention, emotion or belief is reached (a common value for all node states) depends on all these network characteristics. Sometimes silent assumptions are made about the aggregation and timing characteristics. For timing characteristics, often it is silently assumed that the nodes are updated in a synchronous manner, although in application domains this assumption is usually not fulfilled. For aggregation, in social network models usually linear functions are applied, which means that it is often not investigated how a variation of this choice of aggregation would affect the dynamics.

In the modeling and analysis approach used in this paper, a more diverse landscape is covered which is not limited by the fixed conditions on connectivity, aggregation or timing as are so often imposed. For connectivity, both acyclic and cyclic networks are covered here. For aggregation, both networks with linear and nonlinear aggregation are considered and for networks with nonlinear aggregation, networks with logistic aggregation are addressed but also networks with other forms of nonlinear aggregation. Finally, both synchronous and asynchronous timing are provided. The often-occurring use of linear functions for aggregation for social network models may be based on a more general belief that dynamical system models can be analyzed better for linear functions than for nonlinear functions. Although there may be some truth in this if specifically logistic nonlinear functions are compared to linear functions, such a belief is not correct in general. It has been found that also classes of nonlinear functions exist that enable good analysis possibilities when it comes to the emerging dynamics within a network model, thereby among others not using any conditions on the connectivity but instead exploiting for any network its structure of strongly connected components.

In Sec. 6.1 it is shown that for the nonadaptive case this network-oriented modeling approach is equivalent to any dynamical systems modeling approach (Theorem 1 and Corollary 1), and in Sec. 6.2 that for the adaptive case self-modeling networks are equivalent to any adaptive dynamical systems approach (Theorem 2 and Corollary 2). In Sec. 6.3 equilibrium analysis of network models is provided and applied to the model introduced earlier.

6.1. Dynamical Systems and their Canonical Network Representation

Dynamical systems are usually specified in certain mathematical formats; see pp. 241–252 of Ref. 44 for some details. In the first place, a finite set of states (or state variables) $X_1,\dots ,X_n$ is assumed describing how the system changes over time via functions $X_1(t),\dots ,X_n(t)$ of time t. As discussed by Ashby⁴⁴ and Port and Gelder,⁴⁵ a dynamical system is a state-determined system which can be formalized in a numerical manner by a relation (rule of evolution) that expresses how for each time point t the future value of each state $X_i$ at time $t+s$ uniquely depends on s and on $X_1(t),\dots ,X_n(t)$. Therefore, a dynamical system can be described via n functions $F_j(V_1,\dots ,V_n,s)$ for each $X_j$ in the following manner (see also pp. 243–244 of Ref. 44) :

If these functions $F_j$ and the $X_j$ are continuously differentiable (which also implies they are continuous), we call the dynamical system smooth. Suppose such a smooth dynamical system is given. It turns out that it can always be described in a canonical manner by a temporal–causal network model; the argument is as follows. Consider (7) where the functions $F_i$ and $X_i$ are continuously differentiable. In the particular case of s approaching 0 it holds

which due to continuity of the involved functions implies

So, Eq. (7) also holds for $s=0$.

Let $X_j^{\prime }$ denote the derivative of $X_j$ with respect to time. Applying the partial derivative $\partial (..)/\partial s$ on both sides of Eq. (1)

Then it follows Here for the left-hand side, by the chain rule for function composition it holds So, it is found that for all t and s it holds

In particular, this holds for $s=0$; therefore

For a more detailed explanation of the argument for (10), see Appendix  A, where also the differences and relations with Ashby’s approach⁴⁴ are discussed.

Now define the (combination) function $g_j(V_1,\dots ,V_n)$ by

Then it holds

Comparing Eq. (12) to the canonical format (1) in Sec. 3 that defines dynamics of temporal–causal networks, it immediately follows that this matches each other as long as the speed factors $\eta$ and connection weights $\omega$ are set at 1, i.e. from (12) it follows:

with ${\eta }_{Xj}=1$ and $c_{Xj}=g_j$ for all j and ${\omega }_{Xi,Xj}=1$ for all i and j. For an example of this, see Box 1.

This shows that any given smooth dynamical system can be formalized in this canonical manner by a representation in the temporal–causal network format; this notion is described in more detail in Definition 1 and Theorem 1. Note that this also shows theoretically that the use of specific values for speed factors and connection weights is not essential, as they can all be set to 1. However, they still are convenient instruments in the practice of modeling real-world processes.

Definition 1. [Canonical network representation of a smooth dynamical system]. Let any smooth dynamical system be given by

where the functions $F_i$ are continuously differentiable. Then the canonical temporal–causal network representation of it is defined by network characteristics ${\omega }_{Xi,Xj}$, $c_{Xj}$, ${\eta }_{Xj}$ for all i and j with ${\omega }_{Xi,Xj}=1$ for all i and j, $c_{Xj}(V_1,\dots ,V_n)=V_j+{[\partial F_j(V_1,\dots ,V_n,s)/\partial s]}_{s=0}$, ${\eta }_{Xj}=1$ for all j.

This network representation has dynamics induced by the following canonical differential equations for temporal–causal networks

So, by the argument above, the following theorem is obtained:

Theorem 1 (The canonical network representation of a smooth dynamical system). Any smooth dynamical system can be formalized in a canonical manner by a temporal–causal network model called its canonical network representation. Conversely, any temporal–causal network model is a dynamical system model.

As a corollary from Theorem 1 the following well-known result immediately follows.

Corollary 1 (From smooth dynamical system to first-order differential equations). Any smooth dynamical system can be formalized as a system of first-order differential equations.

The latter result was also proven in a different way in pp. 241–252 of Ref. 44. See Appendix  A for some more details.

6.2. Adaptive dynamical systems and their canonical self-modeling network representation

In this section, it is shown how the approach described in Sec. 6.1 can be extended to obtain a transformation of any smooth adaptive dynamical system into a self-modeling network model. Adaptive dynamical systems are usually modeled by two levels of dynamical systems (see Fig. 8).

Here the higher level dynamical system models the dynamics of the parameters $P_{i,j}$ of the lower level dynamical system (the lower level component in Fig. 8) that describes the dynamics of variables $X_i$, for example by

In addition, for the dynamics of the $P_{i,j}$ there will also be a dynamical system (the upper level component in Fig. 8) for $s\ge 0$ :

By applying the argument from Sec. 6.2 to both levels, the following differential equations are obtained covering the entire adaptive dynamical system :

where all $\eta$ and $\omega$ are 1. Recall from Sec. 3 the canonical differential equation (3) that defines a self-modeling network model for the case when self-model states P are introduced for the combination function parameters ${\pi }_Y$ for all base level states Y:

$$\begin{array}{l}dY(t)/dt={\eta }_Y[c_{P_Y(t),Y}({\omega }_{X_1,Y}{X}_1(t),\dots ,\\ {\omega }_{X_k,Y}{X}_k(t))-Y(t)]\end{array}$$

where $P_Y(t)=(P_{1,Y}(t),\dots ,P_{m,Y}(t))$.

The first equation of (15) is (although in a slightly different mathematical notation) equal to Eq. (3), which shows that the former equation defines a self-modeling temporal–causal network model. In this self-modeling network model, the parameters $P_{i,j}$ from the adaptive dynamical system are modeled by (aggregation) self-model P-states within the self-modeling network model for parameters in the combination functions used for states $X_i$ in the base network defined by the states $X_i$. In this way a canonical self-modeling network representation is obtained for the considered smooth adaptive dynamical system; this notion is defined by Definition 2.

Definition 2. [Canonical self-modeling network representation of a smooth adaptive dynamical system]. Let any smooth adaptive dynamical system for $s\ge 0$, $j=1,\dots ,n$, and for $i=1,\dots ,k$ be given by

where the functions $F_j$ and $P_{i,j}$ are continuously differentiable. Then the canonical self-modeling network representation of it is defined by characteristics $\omega$, $\pi$, $\eta$ where all $\omega$ and $\eta$ are 1 and

$$\begin{array}{l}{c}_{Xj}(W_{i,1},\dots ,W_{i,k},V_1,\dots ,V_n)=V_i+{[\partial F_i(W_{i,1},\dots ,W_{i,k},V_1,\dots ,V_n,s)/\partial s]}_{s=0},\\ c_{Pi,j}(W_{i,1},\dots ,W_{i,k},V_1,\dots ,V_n)=W_{i,j}+{[\partial G_i(W_{i,1},\dots ,W_{i,k},V_1,\dots ,V_n,s)/\partial s]}_{s=0}.\end{array}$$

This self-modeling network representation has dynamics induced by the following canonical differential equations:

Thus, by the argument preceeding Definition 2, the following theorem was obtained.

Theorem 2 (The canonical self-modeling network representation of an adaptive dynamical system). Any adaptive smooth dynamical system model can be transformed in a canonical manner into a self-modeling network model called its canonical self-modeling network representation described by Definition 2 above. Conversely, any self-modeling network model is an adaptive dynamical system model. These also apply to higher-order adaptive dynamical systems in relation to higher-order self-modeling networks.

As a corollary it now follows that any adaptive dynamical system can be described by first-order differential equations:

Corollary 2 (From a smooth adaptive dynamical system to first-order differential equations). Any smooth adaptive dynamical system can be formalized as a system of first-order differential equations.

Theorems 1 and 2 demonstrate that a modeling approach based on the self-modeling network format is at least as general as any other adaptive dynamical system modeling approach. Therefore the choice for this format to model adaptive dynamical systems does not introduce any limitation. In particular, it also can be viewed as generalizing the most common types of neural network models.

7. Stationary Point and Equilibrium Analysis for Self-Modeling Networks

In this section, it is shown how stationary point and equilibrium analysis can be performed for self-modeling networks (Sec. 7.1) and how this can be applied to verify the correctness of the implemented self-modeling network model introduced in Sec. 4 compared to its design specifications (Sec. 7.2).

7.1. The general analysis approach

The following types of properties are often considered for equilibrium analysis of dynamical systems in general.

Definition 3 (Stationary point, increasing, decreasing, equilibrium). Let Y be a network state

•	Y has a stationary point at t if $dY(t)/dt=0$,
•	Y is increasing at t if $dY(t)/dt>0$,
•	Y is decreasing at t if $dY(t)/dt<0$,
•	The network model is in equilibrium at t if every state Y of the model has a stationary point at t.

Note that for mathematical analysis of dynamical system models, usually there is an emphasis on equilibrium analysis. However, in many cases no equilibria occur, for example in cases of oscillatory limit cycle behavior. In such cases, still stationary points can be analyzed. In the case, considered in this paper no equilibria occur when the environment is changing all the time.

By considering the canonical network representation of a dynamical system, the above criteria are formulated in terms of the network characteristics: for network models, the following criteria in terms of the network characteristics ${\omega }_{X,Y}$, c_Y, ${\eta }_Y$ can be derived from the generic difference equation (1).³⁸ Let Y be a state and $X_1,\dots ,X_k$ the states connected toward Y. For nonzero speed factors ${\eta }_Y$, the following criteria in terms of network characteristics for connectivity and aggregation apply; here ${aggimpact}_Y(t)=c_Y({\omega }_{X_1,Y}{X}_1(t){\omega }_{X_k,Y}{X}_k(t))$:

•	Y has a stationary point at $t\iff {aggimpact}_Y(t)=Y(t)$
•	Y is increasing at $t\iff {aggimpact}_Y(t)>Y(t)$
•	Y is decreasing at $t\iff {aggimpact}_Y(t)<Y(t)$
•	The network model is in equilibrium at $t\iff {aggimpact}_Y(t)=Y(t)$ for every state Y.

The above criteria for a network being in an equilibrium (assuming nonzero speed factors) depend both on the connections weights ${\omega }_{X,Y}$ used for connectivity and on the combination function c_Y used for aggregation. Note that in a self-modeling network, these criteria can be applied not only to base states but also to self-model states. In the latter case they can be used for equilibrium analysis of learning or adaptation processes.

In particular, a network model with states $X_1,\dots ,X_n$ is in equilibrium if and only if the following n equations (called equilibrium equations) are satisfied:

These equations express relations between values of the states in an equilibrium: they indicate how values in an equilibrium relate to each other and contain as parameters network characteristics ${\omega }_{X_i,X_j}$ and $c_{X_j}$. Sometimes it is possible to solve these equations, for example, when they are linear, or when they are nonlinear Euclidean or geometric equations. When there is no equilibrium, still stationary points for a given state $X_i$ at some time point t can be analyzed based on the above criteria, for example, (local) maxima or minima of the function $X_i(t)$.

The above criteria can be used to verify correctness of (the implementation of) a network model based on inspection of stationary points or equilibria in the following manner.

Verification by checking the criteria through substitution

(1)	Generate a simulation.
(2)	For a sample of states $X_j$ identify stationary points with their time points t and state values $X_j(t)$.
(3)	For each of these stationary points for a state $X_j$ from the chosen sample at time t, identify the values $X_i(t)$ at time t of states $X_i$ among $X_1,\dots ,X_n$ that are connected toward $X_j$.
(4)	Substitute all these values $X_i(t)$ in the criterion ${aggimpact}_{X_j}(t)=X_j(t)$.
(5)	If the equation holds (for example, with absolute deviation $<10^{-2}$), then this test succeeds, otherwise it fails.
(6)	If this test fails, then it should be explored what error is causing this failure and how this error can be corrected.
(7)	If the test succeeds, this contributes to evidence that the implemented network model is correct in comparison with its design specification.

In Sec. 7.2, it is shown in detail how this form of verification by substitution can be applied for an example network.

7.2. Analysis of the introduced neural agent model

The procedure for testing on correctness described in Sec. 7.1 has been applied to the neural agent model introduced in Sec. 4 of this paper. It has been applied for two different scenarios, in each of them for a sample of states covering all levels of the model. One scenario used for this type of test is a scenario where an equilibrium occurs, see Fig. 9. This was achieved by setting the external factors for stimuli and communication enabling to constant values 1 instead of random values as applied in the scenarios in Sec. 5 where no equilibria occur due to this randomness.

The analysis focuses on time $t=1000$ and the chosen sample consists of the following states $X_j$: the second-order self-model state H${}_{WB}$ ($X_{105}$, green line in Fig. 9 ending up around 0.01), first-order self-model state T${}_{\mathrm{\text{exec}}v,B}$ ($X_{95}$, orange line ending up around 0.2), second-order self-model state H${}_{TB}$ ($X_{103}$, red line ending up around 0.75), base states rep${}_{s,B}$ ($X_{31}$, orange line ending up around 0.92), ws${}_{v,B,A}$ ($X_{53}$, purple line ending up around 0.97), intrasyncdet${}_{A,m-v}$ ($X_{19}$, blue line ending up at time 20 around 1), intersyncdet${}_{B,A,m}$ ($X_{21}$, blue line ending up at time 50 around 1), and first-order self-model state W${}_{\mathrm{\text{prep-exec}}v,B}$ ($X_{71}$, blue line ending up at time 500 around 1). To calculate

for each of these chosen states $X_j$, from the network characteristics the states $X_i$ are determined with connections to $X_j$ (see Table 3, second column) and for each of these $X_i$ the weights ${\omega }_{Xi,Xj}(t)$ of the connection from $X_i$ to $X_j$ (see Table 3, middle part), and from the simulation data the values $X_i(t)$ of the $X_i$ at time t (see Table 3, right-hand part).

Based on these connection weights ${\omega }_{Xi,Xj}(t)$ and state values $X_i(t)$ in Table 3, the products ${\omega }_{Xi,Xj}(t)X_i(t)$ are determined (see Table 4, sixth column). Then by applying the combination function from Table 1 on these products, via Eq. (17) aggimpact${}_{Xj}(t)$ is determined for each $X_j$ (see one but last column in Table 4). Finally, this value is compared to the state value $X_j(t)$ of $X_j$ to obtain the $\mathrm{\text{deviation}}={aggimpact}_{Xj}(t)-X_j$ in the last column of Table 4. In Table 4, all absolute deviations are smaller than 10${}^{-6}$. This provides evidence that the implemented model is correct with respect to its design specifications. Note that this provides evidence for correctness of the model in general, not only for this special scenario, as if there were errors, they would most probably also show their effects in this example scenario.

Yet, for still more evidence another scenario has been analyzed as well, where the environmental factors for stimuli and communication enabling do change but in a nonrandom manner; see Fig. 11.

Here, for both agents, stimulus s occurs from time 0 to time 450 and disappears from 450 to 500; this is repeated every 500 time units. Moreover, interaction is not enabled from time 0 to time 50 and is enabled from time 50 to time 400, which is repeated every 400 time units. Due to the changing environment, no equilibrium occurs here. However, there are many cases of (approximate) stationary points. In particular, stationary points have been analyzed as above for a sample of states and time points indicated in the left three columns in Table 5. Here most absolute deviations are $<0.01$. However, there are three of them in the order of 0.02 which is larger than expected for a stationary point. The graph indeed shows that they actually are not approximately stationary points. All in all, also these results provide evidence that the implemented model is correct with respect to its design specifications.

8. Discussion

In this paper, a neural agent model was introduced for the way intrapersonal synchrony and interpersonal synchrony induce behavioral adaptivity between the synchronized persons.^{1,2,3,4,5,6,7,8,9,10} In literature, it was advocated to use a dynamical systems modeling approach to model the complex, cyclical types of dynamics that occur.^11,12 The model presented here is indeed a dynamical system model; moreover, it is multi-adaptive in that the behavioral adaptivity covers both short-term adaptations and long-term adaptations, reflecting short-term affiliation and long-term bonding. The former type of adaptation was modeled using (nonsynaptic) adaptive excitability,^15,16,17,18 whereas for the latter type a more classical synaptic type of adaptation^19,20,21,22 was used. Following the aforementioned literature on synchrony, both types of adaptivity were modeled as driven by the (internally detected) intrapersonal synchrony and interpersonal synchrony for the agent. By also including metaplasticity²³ in the model to control the adaptations in a context-sensitive manner, the agent model became second-order adaptive. The simulations of the model have been performed using the dedicated software environment developed in MATLAB on HP Intel Core i5 and Apple Macbook Pro Intel Core i9 laptops. Execution times per run were less than a minute, for example, on the HP Intel Core i5 with MATLAB 2017a between 40 and 50$$s. The software environment used can be downloaded via URL https://www.researchgate.net/publication/368775720.

Thus, this paper has focused on the emerging and adaptive effects of human social interaction, concerning emerging synchronization and related adaptive affiliation and bonding, with the therapist–client interaction as a central application option. By such simulations, for example, a therapist can get insight about how they can improve the way they interact with their clients and make therapy or counseling more effective. From the scientific perspective, the modeling also contributes formalization to this area of psychology which is almost always addressed in informal manners. The contributed approach can also be used as a solid basis for development of supporting virtual agents in that context.

Synchrony and related patterns in the brain are also analyzed, for example, for atypical brain conditions of subjects such as for PTSD,⁴⁶ epilepsy,^47,48 ADHD,^49,50 or autism.⁵¹ The work presented in this paper distinguishes itself from this in four ways: (1) it abstracts from the specific brain processes and instead focuses on the level of mental processes, (2) it addresses not only the emergence of synchrony but also the causal effects of synchronization on adaptivity of interaction behavior such as affiliation and bonding, (3) it does not address the emergence of synchrony from an objective external observer perspective but from a subjective perspective from the agents themselves, and (4) it focusses on typical instead of atypical conditions of subjects.

We already engaged in computational modeling of synchrony between agents in earlier work.^52,53 However, in the models described there, no (subjective) internal detection of synchrony takes place. Moreover, in the first one no adaptivity was covered,⁵² and in the second one another type of adaptivity was incorporated, namely of internal connections from representation states to preparation states.⁵³ As far as we know, previous work³⁷ describes the only other computational agent model where subjective synchrony detection is addressed. However, by that model no long-term behavioral adaptivity is covered and also no adaptive intrinsic excitability is addressed, whereas both are included in the current model.

Earlier work^1,29 addressing behavioral adaptation due to coordinated actions used a dynamic form of the “bonding based on homophily” principle⁵⁴ was used to model the effect of coordination of emotions and actions on behavioral adaptivity but no (subjective) detection of synchrony was used.

In this paper, also mathematical analysis was addressed for the modeling approach applied. The first type of analysis shows that any smooth adaptive dynamical system has a canonical representation as a self-modeling network. This implies theoretically that the self-modeling network format is widely applicable and that no biases or limitations are introduced by choosing a network modeling approach to design adaptive dynamical system models. In particular, it also generalizes the most common neural system models.

This finding also has been shown in many practical applications varying from biological, cognitive, affective to social processes and their interaction. It is illustrated by many examples in particular in books^38,39 introducing the self-modeling network modeling approach and its applications, a book⁵⁵ focusing on the use of self-modeling network models to handle dynamics, adaptation and control of internal mental models, and a book⁵⁶ focusing on the use of self-modeling network models to model organizational learning processes. Furthermore, stationary point and equilibrium analysis were addressed and applied to the introduced self-modeling network model. These analyses were used as a form of verification of the model, which provided evidence that the implemented model is correct with respect to its design specifications.

Thus, a flexible human-like second-order multi-adaptive neural agent model was obtained for the way in which detected synchrony leads to different types of behavioral adaptivity concerning the short-term affiliation and long-term bonding between the two agents.

A number of aspects that still may be considered relevant are not covered by the model introduced here. One of these aspects is the use of time lags in the process of synchrony detection. This has not been addressed here, but can be a relevant extension of the work reported here. More in general, different combination functions describing methods for synchrony detection may be tried out either with time lags or without time lags. Another relevant aspect that is not addressed in this paper is the role of interruptions or transitions in synchrony and their effect on behavioral adaptivity.

For further work, many more simulation experiments can be designed and conducted, for example to explore the question which types of short-term synchrony are most likely to become translated into long-term benefits for a relationship, or to explore in more detail the roles of intrapersonal synchrony and interpersonal synchrony.

Also, the model can easily be extended to cover interaction between more than two agents. To achieve that, for each additional other agent, within the agent model the number of sensing states can be extended by three of them (for each modality) and similarly three additional representation states and three interpersonal synchrony detection states can be added for each additional other agent. Accordingly, also additional first-and second-order self-model states can be added. This will add complexity to the agent model.

Alternatively, if the model abstracts from the differences between the other agents, then the current agent model without any additional states can be directly applied by using the current sensing, representation, and synchrony detection states and also the self-model states as states aggregating all other agents. This keeps the complexity of the agent model the same, but of course then the model is less context-sensitive as different adaptations to specific other agents are not possible. So, as more often happens, there is a trade-off between the complexity of the agent model and the extent of context-sensitivity here: more context sensitivity comes with more complexity.

Considered from a wider scientific perspective, the model can provide a basis to develop adaptive virtual agents that are able to concentrate on each other by short-term behavioral adaptivity and bond with each other by long-term behavioral adaptivity in a human-like manner. For example, in other work⁵⁷ the focus is on virtual conversational agents and how they can adapt to their human users. In that work, classical learning techniques from AI are used to optimize the agent’s behavior with respect to a given user such as Q-learning, which is not directly inspired or justified by neuroscience.

In contrast, this paper offers approaches such as synaptic plasticity by adaptive connection weights^19,20,21,22 and nonsynaptic plasticity by adaptive excitability thresholds,^15,16 and in addition metaplasticity²³ to control both types of plasticity. As all of these forms of adaptivity are justified in neuroscience literature, this will in principle lead to a more human-like agent model. Nevertheless, it will be interesting to explore in further work how these two different perspectives can benefit from each other.

Concerning the relation of the considered agent model to mechanisms from neuroscience, note that these mechanisms have been incorporated only from an abstracted functional perspective. This means that it cannot be claimed that the model is human-like for the (neuro)physiological level. The latter has been left out of consideration here and would require another research project.

Validation of the model has only been done based on qualitative empirical information from the psychological literature. Dynamic and adaptive patterns have been obtained that are in accordance with that type of empirical information. Due to the lack of quantitative (numerical) empirical information, no quantitative validation has been performed yet. For future research, it is considered to try to acquire such numerical data and then perform quantitative numerical validation by parameter tuning. The way how that can be done is described in Chap. 19 of our 2022 book.⁵⁸

Part of the work presented in this paper was presented in a preliminary form in the AIAI’22 conference and published in their proceedings as a paper⁵⁹ of less than 50% of the length of this paper. This paper is limited to the design of the model and an example simulation. In contrast, the fundamental mathematical analysis of the positioning of the modeling approach based on self-modeling temporal–causal networks in the landscape of adaptive dynamical systems (described in Secs. 6 and 10) is new. It has been shown there that any smooth adaptive dynamical system has a canonical representation as a self-modeling temporal–causal network, which means that the applied modeling approach is universal for smooth adaptive dynamical systems. Moreover, the in-depth verification of the introduced model (described in Sec. 7) is also new. Here, substantial evidence was added that the implemented model is correct with respect to its design specification. Finally, the full specification of the model in Appendix B is new as well.

9. Conclusion

All in all, we achieved the following summarized findings and discoveries:

•	Formalization of the informal domain of human social interaction involving emerging and multiple types of adaptive dynamical system effects is possible.
•	Unifying bridges between causal modeling, network modeling, and dynamical systems modeling are possible.
•	A systematic approach to (higher-order) adaptivity in these different modeling perspectives is possible.

More specifically, the following has been achieved:

•	The notion of canonical temporal–causal network representation for any smooth dynamical system is introduced; it is shown how any smooth dynamical system can be assigned this canonical network representation. This also applies to most common neural network approaches. This creates a bridge between different subdisciplines that are usually kept separate: causal modeling in AI, neural networks in AI, (multidisciplinary) network science, computational science. This enables simulation and analysis of any smooth dynamical system in terms of network concepts.
•	The notion of canonical temporal–causal self-modeling network representation for any smooth adaptive dynamical system is introduced; it is shown how any smooth adaptive dynamical system can be assigned this canonical network representation. This again creates a bridge between different subdisciplines that are usually kept separate and provides a clear way of addressing higher-order adaptivity to any of these: metalevel architectures in AI, causal modeling in AI, neural networks in AI, (multidisciplinary) network science, computational science. This enables simulation and analysis of any smooth higher-order adaptive dynamical system in terms of network concepts.
•	It is shown by mathematical stationary point and equilibrium analysis how verification of correctness of the implemented neural agent model introduced (in comparison to its design specifications) provides more evidence of its correctness.
•	By explaining the full specification of the introduced neural agent model, reproducibility is obtained.