Applying dynamic systems theory and complexity theory methods in psychotherapy research: A systematic literature review

Abstract

Objective

Dynamic systems theory and complexity theory (DST/CT) is a framework explaining how complex systems change and adapt over time. In psychotherapy, DST/CT can be used to understand how a person's mental and emotional state changes during therapy incorporating higher levels of complexity. This study aimed to systematically review the variability of DST/CT methods applied in psychotherapy research.

Methods

A primary studies search was conducted in the EBSCO and Web of Knowledge databases, extracting information about the analyzed DST/CT phenomena, employed mathematical methods to investigate these phenomena, descriptions of specified dynamic models, psychotherapy phenomena, and other information regarding studies with empirical data (e.g., measurement granularity).

Results

After screening 38,216 abstracts and 4,194 full texts, N = 41 studies published from 1990 to 2021 were identified. The employed methods typically included measures of dynamic complexity or chaoticity. Computational and simulation studies most often employed first-order ordinary differential equations and typically focused on describing the time evolution of client-therapist dyadic influences. Eligible studies with empirical data were usually based on case studies and focused on data with high time intensity of within-session dynamics.

Conclusion

This review provides a descriptive synthesis of the current state of the proliferation of DST/CT methods in the psychotherapy research field.

Keywords:

• complex systems
• nonlinearity
• chaos theory
• differential equations
• psychotherapy
• systematic review

Clinical or methodological significance of this article: This review describes the use of mathematical methods in understanding the complexity of psychotherapy dynamics. These methods could be useful tools to build detailed theories about how therapy progresses, including how different aspects change as therapy goes on, which might produce beneficial insight for both therapists and clients. These theories could show for instance how things shift from one phase to another in a particular client's system. However, proving these ideas with real therapy data is tricky due to the limited information that could be gathered during sessions and the unpredictable nature of how the change is measured.

Psychotherapy is a complex activity that includes an ongoing and multifaceted interaction of many variables evolving over time. Both clinical experience and research studies inform us that there are phenomena in the therapeutic process, such as sudden gains (Olthof et al., 2020a), that defy being modeled in a straightforward and linear manner. Change in psychotherapy is often nonlinear (Guastello et al., 2008), and linear models represent rough approximations that cannot do justice to the intricate nature of the therapeutic process. Therefore, several authors have suggested treating psychotherapy as a complex dynamic system that develops over time (e.g., Gelo & Salvatore, 2016; Schiepek et al., 1997).

General systems theory (GST) is an overarching interdisciplinary theory dealing with the study of systems—models where the whole is composed of interconnected parts (Boulding, 1956). GST proposes that systems have common properties across different scientific fields and can be studied generally (von Bertalanffy, 1968). The development of the systems and their components is studied in a subfield dynamic systems theory (DST). DST methods usually describe the development of complex systems in terms of sets of difference (discrete) or differential (continuous) equations. Most suitable differential equations do not have analytical solutions readily available and require the use of advanced mathematical tools to be workable (Strogatz, 2018).

An important group of models derived from GST are the complex systems, the subject of complexity theory (CT), and its related fields (such as chaos theory). While stable systems react to changes in conditions (i.e., perturbation) by returning to a previous state, complex systems behave adaptively, self-organizing into new emerging states (so-called phase transition; Strathern & McGlade, 2014). Chaos theory represents a bridge between deterministic and stochastic conceptions of reality studying unpredictable systems with sensitive dependence on initial conditions which can occur even for deterministic and relatively simple models (Sivakumar, 2004).

The most common paradigm in dealing with complexity is reductionism—attempting to decrease the complexity by studying the parts of the system separately or simplifying the relationships between them (e.g. linearization). While this approach led to several great scientific advances, especially in natural sciences (e.g. classical mechanics in physics), it can lead to a critical loss of information when dealing with more complex systems and their emerging behavior. Evidence of these irreducible complex systems can be found both in natural (Strogatz, 2018) and social sciences (Mellado et al., 2022) creating a need for adoption and focus on complexity theory.

From the DST/CT perspective, psychotherapy can be described via several distinctive but universally recognized complex dynamic system features that include (a) a large number of interrelated components or agents (e.g., therapeutic factors, interacting structures) that can produce emerging self-organizing functioning at the whole-system level; (b) nonlinear temporal dynamics and self-emerging structures while adapting to the environment; and (c) occasionally stochastic nature and an emergent structure that cannot always be traced back to the initial elements of the system (McGill et al., 2020). The eventual outcome of the treatment process may vary greatly even with a slight change in the initial conditions (Schiepek et al., 2017).

In psychotherapy, studying change processes through repeated measurements of variables is essential. This can reveal DST/CT phenomena such as critical fluctuations, sudden gains or losses, and potential chaotic behavior in the therapeutic systems. Macroscopic patterns may emerge, influenced by order and control parameters (the former are measurable quantities characterizing emerging patterns, and the latter are parameters influencing the system’s behavior), which are not visible at the microscopic level. This self-organization process, a core concept of synergetics (Haken & Tschacher, 2017), involves the system evolving towards preferred or stable states represented by attractors. These states are shaped by control and order parameters, guiding the system's dynamics and potentially changing attractors over time (Schiepek et al., 2017). It should also be noted that not all psychotherapeutic change processes must be complex in nature. The “true” behavior may be more regular, linear, and deterministic. However, while some of these systems may be deterministic (i.e., fully predictable from the initial conditions), others may include a stochastic component (i.e., randomness inherent in the system behavior).

In their literature review, Mellado et al. (2022) summarized results obtained from studies employing DST/CT and the DST/CT methods used to operationalize therapeutic change. However, their study (1) was limited to empirical studies only and excluded simulation and mathematical studies (i.e., studies based on simulated data only and studies that only specified a mathematical model), and (2) they did not focus on the specification of the mathematical models used in the primary studies.

Study Aims

A systematic review can bridge the gap between theory and practice, synthesizing literature and providing an overview of DST/CT methods in psychotherapy research. Traditional linear models often fall short of capturing the complexity of psychotherapeutic processes. By reviewing DST/CT methods, we can identify approaches that better capture this complexity within one easily accessible resource for researchers, practitioners, or students, distilling complex mathematical concepts into a cumulative format.

The present systematic review built on Mellado et al.’s (2022) study and aimed to synthesize DST/CT methods employed in existing studies to model empirical or simulated time-series data in the context of psychotherapy. Other descriptive data extracted from each included study were utilized to provide more information about the specific context of the use of these methods. Moreover, the preregistration and online supplemental materials (https://osf.io/k3845/) accompanying this review include detailed and more precise descriptions of methods (e.g., search and selection process, all extracted information) and results (e.g., phenomena, specifications of mathematical methods).

Methods

Selection of Primary Studies

This review adheres to the Preferred Reporting Items for Systematic Reviews and Meta-Analyses (PRISMA) guidelines (Moher et al., 2015). We used the following databases to search for articles in 06/2021: Academic Search Ultimate, Academic Search Complete, PsycINFO, SocINDEX with Full Text, PsycARTICLES, eBook Academic Collection, MathSciNet, and Web of Science. See Supplement 1 for the search string and decisions made during the search process. We considered studies eligible if they (1) employed any DCT/CT modeling approach (only “hard systems” were included, cf. Carey et al., 2015; studies describing qualitative, action-based methodologies termed “soft systems” and studies using only noise reduction or information processing methods were excluded); (2) focused on psychotherapy and counseling phenomena (we excluded studies focusing on psychopharmacotherapy, biological, chemical, and brain research, as well as purely mathematical studies); (3) included either real or simulated data, or at least a precise mathematical model (i.e., we excluded theoretical studies); and (4) focused on adult humans (i.e., we excluded studies focused on children or animals).

To reduce the potential for publication bias, we included gray literature (i.e., unpublished or nonpeer-reviewed studies accessible in EBSCO Open Dissertations). While no restrictions were placed on the date of publication, only articles written in English were included.

The screening process was divided into two phases: abstract and full-text screening. The abstract screening was conducted by the first author. During this phase, only articles that failed to meet the inclusion criteria were excluded; all remaining articles were kept for full-text inspection. The full-text screening was conducted by the first two authors (each screened half of the articles independently; ambiguous articles were discussed by both authors until a consensus was reached).

Data Extraction and Synthesis

We extracted four broad arrays of information from the primary studies: (1) basic descriptive data; (2) analytical/mathematical methods and models used; (3) psychotherapeutic/psychological phenomena that were analyzed; (4) information about empirical data (if collected in the study), and (5) auxiliary information. See Supplement 2 for the detailed list of all aspects that were coded (note that we use a different “numeral system” in the supplementary materials). The extracted information was coded into a spreadsheet file (see Supplement 3). Data extraction was conducted by the first and second authors discussing each included study until a consensus was reached. Extracted content (i.e., topics emerging from the primary studies) was later clustered into more general categories based on content similarity (using a qualitative content analysis)—usually, the resulting categories were highly interwoven and the following division was rather artificial (mainly regarding DST/CT phenomena). During the data synthesis, extracted information was aggregated using a narrative and interpretative synthesis to gain insight across studies for which aggregation was possible and relevant.

Results

Out of the 38,216 articles identified during the search process, we excluded 34,022 articles in the abstract screening phase and an additional 4,040 articles in the full-text screening phase. This resulted in 143 articles that met the initial inclusion criteria. Since this considerably exceeded the number of studies suitable for a narrative and interpretive synthesis, we further decided to exclude all studies that did not include any psychotherapy intervention directly (n = 102), resulting in a final number of 41 studies (see Table I). See Supplement 4 for the flow diagram.

Table I. Descriptive information about all studies included in the review.

ID /#	Reference	Country codes (1b)	Mathematical methods and models (2a); (2b) [(2c)]; (2d), (2e); (2i)	Psychotherapeutic phenomena (3a); (3b) [(3d)]	Software (4b)	Data (5b); (5a); (5c) [(5d)]
1	Baker (2018)	United States of America	1; stable attractors, emotion inertia, critical points, interpersonal influence based on emotional valence [Gottman’s differential equations model]; ODE, first-order; Ch = 0	2; therapeutic relationship as client-therapist affects interaction [cognitive, psychodynamic, emotion-focused therapy—0]	MATLAB, SPSS	1; video-transcripts; 1 [n = 2]
2	Bob (2007)	Czech Republic	3; self-organization of behavior during critical periods, chaotic patterns and state transitions [Lyapunov exponents, false nearest neighbors and embedding dimension]; NA; Ch = 1	1; physiological emotion-related sympathetic activity [hypnosis—0]	Dataplore	1; audio-/video-transcripts, electrodermal aktivity; 0 [n = 2]
3	Brice (2019)	United Kingdom	1; dynamics of treatment pathways [agent-based model, system dynamic model]; ODE, first-order; Ch = 0	1; depression [CBT—1]	R	Only simulated data
4	Buckman et al. (2018)	United States of America	NA; NA [computational physiology model of baroreflex on cardiovascular system]; NA; Ch = 0	1; physiological arousal, baroreflex, heart rate [behavioral biofeedback—0]	NA	1; baroreflex physiological states; NA [n = NA]
5	Burger et al. (2020)	Netherlands	1; rate of change, decay, inertia and influence of other variables [differential equations model, computation model of functional analysis to create a formal theory]; ODE, first-order; Ch = 0	1; anxiety, avoidant behavior, perceived benefits or costs of coping, credibility of catastrophic interpretation [CBT exposure and cognitive reappraisal—0]	R	Only simulated data
6	Fisher et al. (2011)	United States of America	1; degree of order, patterns in dynamic system, rank-order differences, intraindividual variations [spectral analysis]; NA; Ch = 0	1; anxiety [CBT, applied relaxation behavioral therapy—0]	LISREL	2; diary; 0 [n = 33]
7	Fuhriman and Burlingame (1994)	United States of America	NA; embedded patterns of order within and across sessions [sequential analysis of local patterns, correlation integral of global patterns—fractal estimates]; NA; Ch = 1	2; group verbal interactions [general therapy—1]	NA	1; audio-transcripts; 0 [n = NA]
8	Gumz et al. (2014)	Germany	1; topological complexity, degree of synchrony, discontinuous changes [dynamic complexity algorithm—fluctuation + distribution]; NA; Ch = 0.5	3; therapeutic relationship as client-therapist transference-countertransference interactions [psychoanalytic therapy—0]	NA	3; self-report questionnaire; 1 [n = 1]
9	Gumz et al. (2010)	Germany	1; topological complexity, discontinuous changes—abrupt vs. V-shaped changes [dynamic complexity algorithm—fluctuation + distribution]; NA; Ch = 0.5	3; therapeutic relationship as client-therapist transference-countertransference interactions [psychodynamic therapy—0]	NA	3; self-report questionnaire, video-transcripts; 0 [n = 6]
10	Haken and Tschacher (2017)	Germany	1; attractor dynamics, self-organization, agent responses, temporal evolution [synergetic computer algorithm based on differential equations]; ODE, second-order; Ch = 0	1; borderline personality disorder, cognitive-behavioral-experience psychological patterns [general thrapy—0]	NA	Only model for simulation
11	Helmich et al. (2020)	Netherlands	1; sudden changes, relative frequency of other types of therapeutic development change trajectories—null, linear, log-linear, stepped [linear regression on trend]; NA; Ch = 0	1; depression [integrative treatment mostly CBT, psychoeducation, physiotherapy, creative and psychomotor therapy—1]	R	2; self-report questionnaire; 0 [n = 329]
12	Jason et al. (2020)	United States of America	1; temporal evolution, change trajectories [differential equation model—conceptual dynamic influence model]; ODE, first-order; Ch = 0.5	1; alcohol/drug dependence, social support, abstinence self-efficacy, stress [community-based posttreatment support—1]	NA	Only simulated data
13	Kowalik et al. (1997)	Germany	3; critical state transitions, topological complexity, entropy variations in signal, chaoticity, sudden changes, chaos-to-chaos transitions, degree of synchrony [pointwise PD2 algorithm or correlation dimension, entropy rates, Shannon entropy, local Lyapunov exponents, plasticity of chaotic measures—relative-within-session and short-term]; ODE, first-order; Ch = 1	2; therapeutic relationship as client-therapist verbal and nonverbal interactions [solution-oriented therapy—0]	NA	1; video-transcripts; 1 [n = 1]
14	Laramée et al. (2016)	Canada	3; state transitions, fixed cycles, treatment discontinuations [Markov model, microsimulation model]; NA; Ch = 0	1; alcohol consumption, drinking patterns [alcohol dependence pharmacotherapy with psychosocial support—0]	NA	2, 3; daily alcohol consumption; 0 [n = NA]
15	Lauro Grotto et al. (2014)	Italy	3; second-order phase transitions in collective dynamics, degree of social mixing, evolution of the psychological state coupled with the evolution of the group-subjects affinity [agent-based model]; NA; Ch = 1	2; group behavior’s structured patterns based on cognitive-affective evolution [psychoanalytical Bion’s group therapy—1]	NA	Only simulated data
16	Li and Kivlighan (2019)	United States of America	1; mutual influence patterns, multilevel data disaggregation [Gottman’s differential equations model]; ODE, first-order; Ch = 0	2; therapeutic relationship as alliance—client’s perception of session quality [general therapy—0]	Hierarchical Linear Modeling software	3; self-report questionnaire; 1 [n = 180]
17	Liebovitch et al. (2011)	United States of America	1; stable attractors, emotion inertia, critical points, interpersonal influence based on emotional valence [Gottman’s differential equations model]; ODE, first-order; Ch = 0	2; therapeutic relationship as client-therapist affectivity interaction [integrated therapy—0]	MATLAB	Only simulated data
18	Lichtenberg and Knox (1991)	United States of America	2; randomness and disorganization of the group’s interaction [Shannon entropy]; NA; Ch = 0	2; group’s social interaction [Yalom’s group therapy—1]	NA	1; audio-/video-transcripts; 0 [n = 32]
19	Olthof et al. (2020a)	Netherlands	1; phase transitions, system reorganization, degree of destabilization and instability [dynamic complexity algorithm—fluctuation + distribution, peak complexity—relative strength of the highest peak in dynamic complexity profile]; NA; Ch = 0.5	3; depression, client’s cognitive-emotional dynamics [integrative therapy—0]	R	2; self-report questionnaire; 0 [n = 328]
20	Olthof et al. (2020b)	Netherlands	1; sudden changes, critical instability, early warning signals, local highest topological complexity, cumulative complexity peak [dynamic complexity algorithm—fluctuation + distribution, peak complexity—relative strength of the highest peak in dynamic complexity profile]; NA; Ch = 0.5	3; depression, client’s cognitive-emotional dynamics [integrative therapy—0]	R	2; self-report questionnaire; 0 [n = 328]
21	Orsucci et al. (2016)	United Kingdom	3; recurrence patterns, power law, autocorrelation decay [recurrence quantification analysis, autocorrelation analysis, Lempel-Ziv entropy, Shannon entropy, autoregressive models—ARMA, GARCH]; NA; Ch = 0.5	2; physiological client-therapist galvanic/linguistic skin response synchronization [short-term CBT—0]	MATLAB	1; audio-transcripts, galvanic skin response; 1 [n = 1]
22	Paz et al. (2021)	Israel	1; intra- and interpersonal damping—average and change in damping, relation to outcomes [dynamic systems model]; ODE, first-order; Ch = 0	3; therapeutic relationship as client-therapist synchrony in vocal affectivity, emotion self-regulation [short-term psychodynamic therapy—0]	MATLAB	1; audio-transcripts, vocal arousal, self-report questionnaire; 1 [n = 30]
23	Peluso et al. (2012)	United States of America	1; stable attractors, emotion inertia, critical points, interpersonal influence based on emotional valence [Gottman’s differential equations model]; ODE, first-order; Ch = 0	3; therapeutic relationship as client-therapist affectivity interaction [gestalt therapy—0]	MATLAB	Only simulated data
24	Redington and Reidbord (1992)	United States of America	3; phase space reconstruction, flow of trajectories, detection of broad-band noise—power spectrum on Hz, fractal dimension [power spectrum analysis, autocorrelation analysis, fractal dimension]; NA; Ch = 1	1; physiological heart rate, cognitive mental stress [insight-oriented therapy conducted by a psychoanalyst—0]	LabVIEW 2, National Instruments Corp., Austin, TX (signal processing); Systat, Systat, Inc., Evanston, IL; DataDesk, Odesta Corp., Northbrook, IL (vizualization)	1; video-transcripts, heart rate physiology; 0 [n = 1]
25	Reidbord and Redington (1992)	United States of America	3; phase space reconstruction, flow of trajectories, detection of broad-band noise—power spectrum on Hz [power spectrum analysis, autocorrelation analysis]; NA; Ch = 1	1; physiological heart rate, cognitive mental stress [brief insight-oriented therapy conducted by a psychoanalyst—0]	SYST AT 5 (Systat); StatView II (Abacus Concepts); Odesta's DataDesk; LabVIEW 2; Mathematica (Wolfram Research); MathCad (MathSoft); Abacus Software's MacSpin	1; video-transcripts, heart rate physiology; 0 [n = 1]
26	Sabelli et al. (1990)	United States of America	3; state transitions—bipolarity [bifurcation model]; NA; Ch = 1	1; bipolar disorder, emotional states [psychotherapy + pharmacotherapy—0]	NA	2; self-report questionnaire; 0 [n = 3]
27	Schiepek et al. (2016a)	AUT/Germany	1; state dynamics, order transitions [dynamic complexity algorithm—fluctuation + distribution]; NA; Ch = 0.5	3; resilience or coping, autonomy, withdrawal, stress [integrated mainly systemic therapy—0]	Synergetic Navigation System (integrated analyses)	2; self-report questionnaire; 0 [n = 1]
28	Schiepek et al. (2014)	Austria	1; complexity, self-organization, local peaks of complexity, critical fluctuations, order transitions [dynamic complexity algorithm—fluctuation + distribution]; NA; Ch = 0.5	3; obsessive compulsive disorder, client’s cognitive-emotional dynamics [CBT—1]	Synergetic Navigation System (integrated analyses)	2, 3; self-report questionnaire; 0 [n = 23]
29	Schiepek et al. (2017)	Austria/Germany	3; synergetics core concepts—self-organized thresholds, synchronization patterns, phase transitions, chaotic dynamics [Schiepek’s differential equations model]; ODE, first-order; Ch = 1	3; client’s cognitive-emotional dynamics [general therapy—0]	MATLAB	Only simulated data
30	Schiepek et al. (2016b)	Austria/Germany	3; synergetics core concepts—self-organized thresholds, synchronization patterns, phase transitions, chaotic dynamics, dynamic patterns with and without the Verhulst-driven eigendynamics of the Emotions order parameter [Schiepek’s differential equations model, Verhulst-driven eigendynamics]; ODE, first-order; Ch = 1	3; client’s cognitive-emotional dynamics [general therapy—0]	MATLAB	2; self-report questionnaire; 0 [n = 1]
31	Schiepek et al. (1997)	Germany	3; topological complexity, chaoticity and sensitivity to initial conditions, nonperiodicity and nonpredictability of time-signal—correlation decay—detection of broad-band noise—power spectrum [sequential plan analysis, Fast Fourier Transformation on power spectra, autocorrelation analysis, determinism test, correlation dimension, largest Lyapunov exponents, local largest Lyapunov exponents]; NA; Ch = 1	2; therapeutic relationship as a client-therapist interaction of cognitive-emotional schemata—verbally and nonverbally communicated self-presentations [solution-oriented brief therapy—0]	NA	1; video-transcripts; 1 [n = 1]
32	Schiepek and Strunk (2010)	Austria	1; critical instability and fluctuations, phase transitions [dynamic complexity algorithm—fluctuation + distribution, permutation entropy, Lyapunov exponent]; NA; Ch = 1	3; obsessive compulsive disorder, client’s cognitive-emotional dynamics [general therapy—0]	NA	2, 3; self-report questionnaire; 0 [n1 = 94, n2 = 1]
33	Schiepek et al. (2013)	Austria	1; pattern complexity of signal, fluctuation —amplitude and frequency of changes in a signal— and distribution —scattering of values or system states within range of possible states, order transitions, critical instability [dynamic complexity algorithm—fluctuation + distribution]; NA; Ch = 0.5	3; obsessive compulsive disorder, client’s cognitive-emotional dynamics [CBT—0]	R, SPSS	2; self-report questionnaire; 0 [n = 18]
34	Schöller et al. (2018)	Austria/Germany	3; synergetics core concepts—self-organized thresholds, synchronization patterns, phase transitions, chaotic dynamics, critical instability [Schiepek’s differential equations model extended to control parameter dynamics]; ODE, first-order; Ch = 1	3; client’s cognitive-emotional dynamics, personality development [general therapy—0]	psysim.at web software	2; self-report questionnaire; 0 [n = 941]
35	Sharman (1998)	Canada	1; attractor states —repetitive defensive personality patterns—, repellor states —dissipating patterns—, state transitions and phase shifts, density of recurrence, anxiety as a control parameter [autocorrelation lag analysis, dynamic systems state-space model]; NA; Ch = 0	1; cognitive-affective personality components interactions controlled by anxiety levels [short-term psychodynamic therapy—0]	MacSHAPA program	1; video-transcripts; 0 [n = 1]
36	Stapelberg (2016)	Australia	2; spectral power density [Fast Fourier Transformation, autoregression analysis, detrended fluctuation analysis]; NA; Ch = 0	1; physiological heart rate variability, depression, anxiety [CBT—0]	MATLAB, SPSS	1; heart rate physiology, self-report questionnaire data; 0 [n = 88]
37	Tschacher and Haken (2019)	Switzerland	3; attractors, state evolution, entropy, potential function, degree of synchrony [Fokker-Planck equation]; PDE, second-order; Ch = 0	2; physiological nonverbal synchrony [general therapy—0]	MEA (motion energy analysis software)	1; heart rate physiology; 1 [n = 2]
38	Tschacher et al. (2015a)	Switzerland	3; order, invariant lagged relationships, phase space, attractor landscape, relaxation time to resettle after perturbation, damping influence [Landberger’s order, vector autoregression, Fokker-Planck equation]; PDE, second-order; Ch = 0	2; therapeutic relationship as client-therapist experiential, behavioral, and structural alliance [schema-centered therapy—0]	MEA (motion energy analysis software)	Only simulated data
39	Tschacher et al. (2015b)	Switzerland	3; pattern formation, phase transitions, stability and change [Fokker-Planck equation]; PDE, second-order; Ch = 0	2; therapeutic relationship as client-therapist experiential alliance—communion and agency –, or as behavioral alliance—nonverbal synchrony [person-centered Rogerian therapy—0]	MEA (motion energy analysis software)	1; video motion-energy-data; 1 [n = 1]
40	Tschacher et al. (1998)	Switzerland	1; order, self-organization, pattern formation [Landsberger’s order analysis, Shannon entropy]; NA; Ch = 0	3; therapeutic relationship as client-therapist alliance—progress within and outside of therapy setting [CBT, insight-oriented, client-centered, heuristic, schema therapy—0]	NA	3; self-report questionnaire; 0 [n = 28]
41	Walter et al. (2010)	Germany	1; topological complexity, unsteady development of changing complexity [dynamic complexity algorithm—fluctuation + distribution]; NA; Ch = 0.5	2; emotional abstraction patterns—relaxing, reflecting, experiencing, connecting [short-term solution-oriented therapy—0]	NA	1; video-transcripts; 0 [n = 1]

Note. (1b): country codes according to the affiliation of the first author; (2a): 1 = deterministic methods, 2 = stochastic methods, 3 = hybrid methods (a combination of deterministic and stochastic methods); (2b): list of DST/CT phenomena; (2c): list of DST/CT methods; (2d): ODE = ordinary differential equations, PDE = partial differential equations; (2e): order of differential equations; (2i): 0 = chaos not present, 0.5 = chaos theoretically, 1 = chaos present computationally; Hz = Hertz; (3a): 1 = the main phenomena are therapeutic outcomes (or semioutcomes in time-series sense)—focusing on the effect of therapy; 2 = the main phenomena are therapeutic processes—focusing on how the therapy develops regardless of the effectiveness of successful therapy; 3 = a combination of both process and outcome focuses; (3b): list of primary phenomena; (3d): list of psychotherapeutic intervention types—distinction between group = 1 and individual = 0 therapy; CBT = cognitive behavioral therapy; (4b): used software; (5b): 1 = micro (processes during therapy session, measured in seconds/minutes); 2 = meso (processes during daily life of patients, measured in days); 3 = macro (processes during therapy session, measured in weeks/sessions); (5a): list of empirical data types; (5c): 0 = individual data, 1 = dyadic data; (5d): clients’ sample size; NA = not available information.

Descriptive Information

The included studies were published between 1990 and 2021 with a median of 2015. Interestingly, the publication dates have a bimodal distribution, as no studies were published between 2000 and 2005. All studies focused on the application of models to time-series data. Out of the total sample of 41 studies, 24 used their own empirical data, nine studies used only an empirical example from a previous study, and eight studies were purely in silico studies. One study used simulations without properly specifying the mathematical model (Laramée et al., 2016, #14). Six of the empirical studies also included simulations. Two studies were neither empirical nor a simulation study but were focused on the specification of the mathematical models only (Haken & Tschacher, 2017, #10; Tschacher et al., 2015a, #38). Twenty-one of the empirical studies follow the procedure of specifying a research question first, then collecting data, and finally selecting a suitable DST/CT method to answer the research question. However, the remaining 12 studies specified a theoretical and mathematical model first and then used empirical data to showcase the model. Moreover, out of the 21 “data first, model after” studies, 7 used DST/CT tools only to create variables for further linear, noncomplex analyses, such as means-comparison (Helmich et al., 2020, #11; Lichtenberg & Knox, 1991, #18; Olthof et al., 2020a, #19, 2020b, #20; Schiepek et al., 2014, #28; Stapelberg, 2016, #36; Tschacher et al., 1998, #40).

Furthermore, 21 studies contained comprehensive theoretical and conceptual explanations of complexity-related concepts or methods used, while the remaining 20 were written in a more applied manner already assuming certain prior knowledge of DST/CT by the reader. See Table I and Supplement 5 for more detailed descriptive information about all included studies (specifically see Supplement 5a for a distribution of studies divided according to countries of origin and year of publication and Supplement 5b for a figure showing in which journals were the included studies published).

Analytical/Mathematical Methods and Models

Analyzed complexity phenomena

Temporal dynamics are a necessary prerequisite for any dynamic system (Ryo et al., 2019). Several authors described the temporal evolution of the system variables in the included studies and/or focused on specific trajectories of the system (n = 13) or focused on the decrease in the amplitude of the time signal oscillation (n = 4) due to the decay of change over time or due to external influences, such as intrapersonal damping.

Self-organization is a topic in complex dynamic systems and represents the spontaneous emergence of organized behavior in response to certain conditions (Heylighen, 2008). Eight studies directly targeted self-organization, five studies focused on the identification of stable attractors (i.e., a tendency of the system to remain within a certain range of values even after perturbation), and four studies focused on the identification of the emergence of order from disorder unfolding in time (e.g., embedded patterns of order within and across psychotherapeutic sessions as a counterpart to complexity; Heylighen, 2008).

Emergence is a term tightly interconnected with self-organization (being a subset of it). In this review, we divided emergence into several categories: the emergence of phases with possible phase transitions (n = 12), the emergence of critical instability and sudden shifts, such as detecting early warning signals (n = 11), the emergence of synchronous states, such as the degree of synchrony between the client and the therapist (n = 7), and the emergence of chaotic patterns with sensitivity to initial conditions (n = 6). For more detailed information about the extracted DST/CT phenomena and the number of studies that used them, see Supplement 6 (text) and for tabulated data specifically Supplement 6a.

Mathematical methods

Of the 41 studies, 22 employed purely deterministic DST/CT methods, 14 combined deterministic approaches with some stochastic elements, and two studies used only stochastic methods (Lichtenberg & Knox, 1991, #18; Stapelberg, 2016, #36).

As far as the particular techniques used to analyze the DST/CT phenomena are considered, these included agent-based models (n = 2, see Supplement 7a), methods for modeling stochastic processes (i.e., spectral analysis, n = 4; Shannon entropy, n = 4; autocorrelations, n = 5, see Supplement 7b), methods measuring dynamic complexity (n = 10, see Supplement 7c), and measuring chaoticity (n = 4, see Supplement 7d). Agent-based models were used to simulate interactions between individual agents in a given system over time (Brice, 2019, #3; Lauro Grotto et al., 2014, #15). Methods dealing with stochastic processes were used to describe the evolution of a system over time and identify patterns in the data. Dynamic complexity studies time-dependent fluctuations and changes in distributions and was used as a pragmatic alternative to other measures of time-series complexity that require larger numbers of data points (e.g., correlation dimension, Schiepek et al., 2017, #29). The chaoticity and stability of the system are usually measured by the local largest Lyapunov exponent. These methods are often used to transform raw data into time series of complexity/instability measures, enabling better identification of phase-space/order transitions. Typically, studies used a combination of methods rather than a single method. For more detailed information about the methods, see Supplement 7 in general.

Mathematical models based on differential equations

Sets of differential equations are used to describe time-dependent complex system behavior. A differential equation describes a relationship between dependent and independent (often time) variables and contains derivatives of dependent variables (Zill, 2012). Numerical methods allow us to conduct simulated experiments even without collecting empirical data. Solving a differential equation requires knowledge of the initial (value of a system at T = 0 at the initial state before time-evolution) and boundary conditions (if known, represent constraints for the solution of the equation).

In the present review, 16 studies specified a model based on differential equations. Of these, 13 studies used ordinary differential equations (ODEs; comprising unknown functions of a single independent variable, unknown function derivatives, and other functions of the independent variable), while three studies used partial differential equations (PDEs; comprising unknown functions of more than one independent variable and their partial derivatives to model multidimensional systems) (Tschacher & Haken, 2019, #37; Tschacher et al., 2015a, #38, 2015b, #39). Except for one study, all ODEs were first-order differential equations. All PDEs were second-order differential equations. Numerical methods are often applied to solve ODEs and PDEs because their analytical solution is usually unavailable. For example, Runge‒Kutta methods (Butcher, 1996) are used to estimate a numerical solution to ODEs by applying discretization of the temporal derivatives. None of the included studies in the present review provided an analytical solution for their set of equations.

The three most influential models employing differential equations included (1) the adjusted Gottman et al.’s (2002) model of dyadic interaction (Liebovitch et al., 2011, #17) leading to stable attractors (Schiepek et al., 2016a, #27), see Supplement 8a; (2) Schiepek’s model of psychotherapy including problem intensity, success in therapy, motivation to change, emotions, and insight (see Supplement 8b); and (3) the Fokker-Planck equation (Tschacher et al., 2015b, #39), see Supplement 8c.

After 2019, more models outside of these three aforementioned branches of models emerged (see Supplement 8d for details). Brice (2019, #3) focused on a system dynamics model with inputs from an agent-based model to describe the dynamic evolution of a system of depressed patients flowing from general practitioners’ care to visits of psychotherapists and psychiatrists and back in several feedback loops. Burger et al. (2020, #5) created a computational model of eight differential equations of one hypothetical patient to create a formal idiographic theory, modeling the potential effects of two possible interventions on the vicious circle of discriminant stimuli, catastrophizing, panic, avoidance, and credibility of catastrophic interpretations. Jason et al. (2020, #12) focused on describing dynamic interrelations between social support, abstinence self-efficacy, and stress using three differential equations. Paz et al. (2021, #22) employed a combination of differential equations with mixed models to simultaneously model the first-order change in attractor patterns while allowing variations at multiple levels, which is important when using data with a nested structure. This approach can be further utilized, for instance, to directly predict therapeutic outcomes.

Emergence of chaotic patterns

Based on the included studies, psychotherapy seems to be a discipline where the emergence of chaos is possible. Even though dynamic systems modeled in the included studies were often stable and predictable (e.g., leading to fixed-point attractors), several studies investigated the emergence of chaotic patterns (n = 14). The emergence of deterministic chaos mainly in the unpredictable client-therapist interaction (Kowalik et al., 1997, #13; Schiepek et al., 1997, #31, 2016b, #30) was inferred from the positive Lyapunov exponents (Bob, 2007, #2; Schiepek & Strunk, 2010, #32), bifurcations creating chaotic states (Redington & Reidbord, 1992, #24; Sabelli et al., 1990, #26; Schiepek & Strunk, 2010, #32), visual inspection of phase-space portraits (Reidbord & Redington, 1992, #25), phase transitions, multistability, chaos-to-chaos phase transitions (Schiepek et al., 2017, #29; Schöller et al., 2018, #34), and the emergence of nonlinear patterns of global complexity with a low-dimensional chaos-like system or fractal dimension (Fuhriman & Burlingame, 1994, #7; Lauro Grotto et al., 2014, #15). Furthermore, 10 other studies suggested that chaos was potentially present but did not test it directly.

Visualization

The included studies used various methods to visualize the results. Visualizations were often provided for both raw data and for DST/CT phenomenon measurement trajectories. Types of visualizations used in the included studies comprised line graphs (e.g., n = 17 raw data, or n = 10 DST/CT data), phase-space portraits (n = 13), heatmaps and/or resonance diagrams (n = 7), cluster dendrograms and/or bifurcation diagrams (n = 4), recurrence plots (n = 3), and conceptual models. See Supplement 9 for more information about the visualization methods used.

Psychotherapeutic/Psychological Phenomena and Interventions

The studies focused on therapeutic outcome (n = 14), therapeutic process (n = 13), or a combination of process and outcome (n = 14). In summary, the psychotherapeutic phenomena suitable for DST/CT methods that were investigated in the included studies could be clustered into more general categories of (1) psychotherapeutic relationship (most investigated, e.g., alliance, dyadic reactions, transference), (2) group processes (e.g., member interactions), (3) physiological processes (e.g., heart rate variability), (4) therapeutic destabilization of pathological states, and (5) specific syndromes (e.g., depression, bipolar disorder) and ideographically approached phenomena (e.g., client cognitive–affective dynamics). A detailed description of extracted psychotherapeutic/psychological phenomena and interventions is provided in Supplements 10a and 10b.

Studies that Used Empirical Data

In this section, we focus in more detail on studies that used empirical data (i.e., we omit the nine simulation studies). Overall, nine studies used only an empirical example based on data from another study (Buckman et al., 2018, #4; Fisher et al., 2011, #6; Fuhriman & Burlingame, 1994, #7; Laramée et al., 2016, #14; Sabelli et al., 1990, #26; Schiepek & Strunk, 2010, #32; Schiepek et al., 2016a, #27; Schöller et al., 2018, #34; Tschacher et al., 2015b, #39) and not a proper cross-validation. In terms of the type of data and data granularity considered, we extracted the smallest unit of measurement for each variable in the included studies and divided them into micro, meso, and macro levels. At the micro level, 16 studies analyzed data with high time intensity focusing on the processes during therapeutic sessions (e.g., heart rate variability and verbal utterances). At the meso level, 9 studies analyzed data focusing on processes in patients’ daily life, typically using data from daily administered self-report questionnaires. At the macro level, four studies analyzed data focusing on the overall trajectories of the therapeutic change based on session-by-session or weekly measurements. Additionally, three studies combined meso and macro levels (e.g., using two different empirical examples with different data granularity). Overall, nine studies used dyadic data, while the remaining studies used individual data (typically based on patients’ perspectives). See Supplement 11 for more information.

Auxiliary Information

This review identified 23 studies that claimed to bring novelty to the field. The list of innovations demonstrates that the field of DST/CT application in psychotherapy research is still at its beginning—nearly half of the studies make novel contributions (see Supplement 12a for information about claimed novelty of the studies). The most often used software was MATLAB. However, R seems to be increasing in popularity in recent years and is slowly replacing the usage of MATLAB in this field. Open data or code are not truly common thus far in terms of open and reproducibility science. For more detailed information about software used and transparency see Supplement 12b. We also gathered several claimed limitations from the included studies, see Supplement 12c.

Discussion

This review narratively synthesized descriptive information from 41 psychotherapy-related studies applying methods from dynamic systems theory and complexity theory (DST/CT) frameworks within simulated or real time-series data. The present review aimed to provide information about the current state of adoption of these methods and their application in the field. Our review shows that these techniques have provided powerful tools for analyzing psychotherapy as a complex dynamic system and understanding the emergent behavior of these systems.

This review showed that the topic of dynamic modeling and complex systems in psychotherapy is currently understudied. We observed that the range of dynamic modeling methods used before 2019 was rather narrow, as only three models were developed. They utilized only one of the following specific models or their variations: (1) describing the therapeutic relationship as a set of influence functions between client-therapist dyad (based on an adapted marital dyads model, Gottman et al., 2002), (2) coevolution of key therapeutic variables and their control parameters, and (3) as a synchrony.

Recently, several other models and applications have emerged, driven by readily accessible computational resources and advanced easy-to-use modeling software. The most recent of the included studies focused mainly on combining multiple methods (Brice, 2019, #3; Paz et al., 2021, #22), constructing formal theories (Burger et al., 2020, #5), continuing with formalizing specific dynamics of selected variables inspired by the tradition of Liebovitch et al.’s (2011, #17) approach and specifying their own model (Jason et al., 2020, #12). Other models, such as the model of couples therapy adding another layer to client-therapist interactions, remain to be operationalized and tested in the context of psychotherapy.

New researchers are encouraged across other countries to continue in this regard. Researchers collecting specific data or wanting to implement complexity in their analyses might want to use the shortcut provided by our review to quickly learn about what has been done, which research questions have been pursued, and what methods may be feasible for their specific data (see Supplement 3). Further studies should examine the use of other DST/CT methods, such as discrete event simulation, optimization algorithms, evolutionary algorithms, cellular automata, and swarm intelligence, within the context of psychotherapeutic interventions. For a scoping review of the use of machine learning in psychotherapy research, see Aafjes-van Doorn et al. (2021).

The importance of computational accessibility is demonstrated by the fact that the majority of computational studies were performed on personal computers and not in supercomputer centers. The importance of the easy accessibility of the simulation capability is demonstrated by the frequent use of high-level languages such as MATLAB or R for simulations as opposed to low-level C or Fortran. Additionally, the authors of the included studies usually did not explicitly mention specific numeric methods to solve the equations, such as linear multistep methods or Runge‒Kutta methods (Butcher, 1996), but used readily available tools, e.g., MATLAB.

Authors usually simulated their data using differential equations, plotted the attractor dynamics, or computed characteristics of the dynamical systems, such as the largest Lyapunov exponents or fractal dimensions. However, even though psychotherapy researchers already assumed self-organization and the emergence of complex behavior of psychotherapeutic systems, many of the identified studies remained in the computational and theoretical framework without cross-validation using empirical data. A similar situation has been observed in the field of health research (Rusoja et al., 2018). Validating the assumptions, as well as tentative conclusions of the existing studies, using real-life data thus remains another challenge in the DST/CT field.

Additionally, a relatively large number of studies (approximately 50%) employed DST/CT methods only to provide transformation of raw data into complexity measures, further imputed in more linearized analytical methods to compare means or seek relationships between specific variables, usually in a multilevel mixed model framework. This is supported by the dynamic models in the individual studies where estimated attractors were often stable, leading to fixed point or limit cycle equilibriums. However, linearization is a generally used practice when dealing with nonlinear systems (cf. Cheng et al., 2010), the full potential of DST/CT methods thus often remains unemployed in a psychotherapeutic context. When used as a full-bodied method, DST/CT may offer an apparatus to create formal psychotherapeutic theories describing the evolution of state variables given the evolution of control parameters enabling transitions between qualitatively different system phases. Importantly, several studies have identified chaotic behavior in psychotherapy data with an important example of Kowalik et al. (1997, #13), who also showcased chaos-to-chaos transitions. Additionally, while deterministic chaos has gained attention in the past, its prominence as a serious topic of study has waned in recent years, with a few exceptions. Empirical validation of chaotic dynamics in actual therapeutic time series remains limited mainly because of limited data resolution (i.e., number of measurement points) and the stochastic nature of measurement issues. However, future advancements in data collection methods may provide opportunities for a more nuanced exploration of chaotic dynamics in the therapeutic processes. Moreover, theoretical explorations and simulations of deterministic chaos in psychotherapy have already yielded valuable insights.

The focus on case studies and time-intensive collection of within-session data, both physiological and verbal, was expected given that the applied DST/CT methods usually require large numbers of measurement points. However, some methods, such as dynamic complexity measures that were already ideographically implemented in clinical practice, were adopted even for daily-based questionnaire data using dynamic complexity instead of raw time series to inform about sudden changes and critical instability (Olthof et al., 2020a). The DST/CT methods, therefore, also demonstrated the possibility of approaching psychotherapeutic phenomena from the ideographic perspective while remaining in the quantitative framework. In such cases, they allowed for the computation of simple complexity estimations based on fluctuation and distribution of time series even of the single item in as few as 20 observations (Schiepek & Strunk, 2010). Moreover, dynamic complexity was validated as a measure of critical instability and phase transitions using correlations with other similar measures, such as permutation entropy or Lyapunov exponents. This is a pragmatic tool accessible even when naturalistically less rich data are analyzed and could be easily employed by other psychotherapy researchers.

This review aimed to provide a description of the applied mathematical DST/CT methods and summarize specific models used in the studies. The review suggests that the application of DST/CT methods in psychotherapy research is still in its infancy, with many studies demonstrating novel contributions but seldom cross-validating the findings except by using empirical examples. Moreover, if we cluster the authors’ affiliations to institutions on the country level, two major clusters emerge: USA and German-speaking countries (i.e., Germany, Austria, and Switzerland). Furthermore, it seems that the majority of included studies are influenced by a small number of influential researchers who coauthored many of these studies. For instance, one of the researchers coauthored 14 of the included studies. However, the validation and calibration of these methods together with seeking good-practice guidelines is an important next step in the development of the field.

The present review has several limitations (see Supplement 13 for more details). The review does not include SCOPUS as a source of studies (because of an excessive number of search results) and the very recent developments in the field (e.g., dynamic latent-class structural equation modeling or Markov switching models, Flückiger et al., 2022) because the search was conducted in 2021 and because several relevant keywords might have been excluded for producing a large number of false positives (e.g., “synchrony”) or omitted unintentionally. For specific reviews targeting “synchrony in psychotherapy,” see Kleinbub (2017) or Koole and Tschacher (2016). Moreover, studies targeting a specific subset of complexity theory, e.g., network modeling, or studies targeting other than human adult populations were omitted arbitrarily, limiting the potential prevalence of other eligible models. Finally, the abstract screening was conducted only by a single reviewer, the methodological quality of the included studies was not assessed, and some variables were only subjectively evaluated (e.g., theoretical comprehensiveness of included studies). Future reviews should focus on novel models and approaches such as computational neuroscience.

In conclusion, to utilize the descriptive nature of the present review, characteristics of individual studies were dummy coded according to several criteria to provide a quick combination of identifiers that the readers might want to filter for when seeking the exact inspiration for their DST/CT research question. A relatively large number of studies provided not only a set of mathematical formulae but also comprehensive theoretical or conceptual explanations. This might clarify the fundamental principles of DST/CT and the intended usage in psychotherapy.

CRediT Authorship Statement

(1) Adam Klocek: funding acquisition, project administration, conceptualization, data curation, formal analysis, investigation, methodology, writing—original draft, writing—review & editing; (2) Jan Premus: conceptualization, formal analysis, methodology, supervision, writing—review & editing; (3) Tomáš Řiháček: conceptualization, supervision, writing—review & editing.

Disclosure Statement

No potential conflict of interest was reported by the author(s).

Supplemental Data

Supplemental data for this article can be accessed online at https://doi.org/10.1080/10503307.2023.2252169.

Additional information

Funding

This study was supported by Specific University Research Grant No. MUNI/A/1187/2021 provided by the Czech Ministry of Education, Youth, and Sports.