电场作用下神经元网络的奇异态

Chimera states in a neuronal network under the action of an electric field

阅读

• 下载地址：https://chinazhang-my.sharepoint.com/:b:/g/personal/rep_rebirth_zh_zhangwenhao_icu/EXsFqKmixQ9JhU9sPBRD0PYBxMFC3B3SfkBv3fVxRFwA0w?e=QP9VNK

module HR
    implicit none
    real,parameter :: h=0.01
    integer,parameter :: N=100,MaxT=170000
    integer,allocatable :: neighbour_matrix_2(:,:),neighbour_matrix_2P(:,:)
    real :: x(N),y(N),z(N)
contains

    subroutine x0_y0_z0()
        implicit none
        integer :: k
        real :: x1,x2,x3
        call random_seed()
        do k=1,N,1
            call random_number(x1)
            call random_number(x2)
            call random_number(x3)
            x(k)=0.001*(k-N/2.0)+x1
            y(k)=0.002*(k-N/2.0)+x2
            z(k)=0.003*(k-N/2.0)+x3
            write(10,*) x(k),y(k),z(k)
        end do
    end subroutine x0_y0_z0

    subroutine fnf(p)
        implicit none
        integer :: i,j,p,t
        real :: func_x(N),func_y(N),func_z(N),func_e(N),a,b,d,k1,k2,k3,k4,II,JJ,C,r,s,x_i0,E_ext,x_s,&
                &lambda,thet_s,C1,C2,E(N)
        k1=0.0
        k2=0.0
        k3=0.0
        k4=0.0
        a=1.0
        b=3.0
        d=5.0
        r=0.001
        s=4.0
        x_s=2.0
        x_i0=-1.6
        lambda=10.0
        thet_s=-0.25
        E_ext=0.0
        E=0.0
        II=4.0
        do t=1,MaxT,1
            if(mod(t,5000)==0) then
                write(*,*) t
            end if
            do i=1,N,1
                JJ=0.0
                C1=0.0
                C2=0.0
                C=0.0
                do j=1,N,1
                    JJ=JJ+neighbour_matrix_2(i,j)*(x(j)-x(i))
                    C1=C1+neighbour_matrix_2P(i,j)*(1.0/(1.0+exp(-lambda*(x(j)-thet_s))))
                    C2=C2+neighbour_matrix_2(i,j)*(1.0/(1.0+exp(-lambda*(x(j)-thet_s))))
                end do
                C=C1-C2
                func_x(i)=x(i)+h*(y(i)-a*x(i)**3+b*x(i)**2-z(i)+II+k3*JJ+(k4/(2.0*p-2.0)*(x_s-x(i))*C))
                func_y(i)=y(i)+h*(1.0-d*x(i)**2-y(i)-k1*E(i))
                func_z(i)=z(i)+h*(r*(s*(x(i)-x_i0)-z(i)))
                !func_e(i)=E(i)+h*(k2*y(i)+E_ext)
                if(mod(t,50)==0) then
                    write(20,*) i,t,x(i)
                end if
                if(i==5.and.t>70000) then
                    write(30,*) t,x(i)
                    write(70,*) x(i),y(i),z(i)
                    write(*,*) t,x(i)
                end if
                if(i==50.and.t>70000) then
                    write(40,*) t,x(i)
                end if
            end do
            x=func_x
            y=func_y
            z=func_z
        end do
        return
    end subroutine fnf


    !建立网状结构

    !k=2
    subroutine neighbour_2(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix_2(N,N))
        !初始化网络矩阵
        neighbour_matrix_2=0
        do i=1,N-1,1
            do j=i+1,N,1
                if(abs(i-j)<=K/2.or.abs(i-j)>=N-K/2) then
                    neighbour_matrix_2(i,j)=1
                    neighbour_matrix_2(j,i)=1
                end if
            end do
        end do
        !输出矩阵
        do i=1,N,1
            do j=1,N,1
                write(50,"(I2)",advance='no') neighbour_matrix_2(i,j)
            end do
            write(50,*)
        end do
        return
    end subroutine neighbour_2
    !k=2*p+2
    subroutine neighbour_2P(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix_2P(N,N))
        !初始化网络矩阵
        neighbour_matrix_2P=0
        do i=1,N-1,1
            do j=i+1,N,1
                if((abs(i-j)<=K/2.or.abs(i-j)>=N-K/2).and.abs(i-j)/=1) then
                    neighbour_matrix_2P(i,j)=1
                    neighbour_matrix_2P(j,i)=1
                end if
            end do
        end do
        !输出矩阵
        do i=1,N,1
            do j=1,N,1
                write(60,"(I2)",advance='no') neighbour_matrix_2P(i,j)
            end do
            write(60,*)
        end do
        return
    end subroutine neighbour_2P

end module HR



program main
    use HR
    implicit none
    integer :: p
    p=40
    open(10,file="data_x0_y0_z0.txt")
    open(20,file="data_i_t_x.txt")
    open(30,file="data_t_x5.txt")
    open(40,file="data_t_x50.txt")
    open(50,file="neighbour_matrix_2.txt")
    open(60,file="neighbour_matrix_2P.txt")
    open(70,file="x.txt")
    call x0_y0_z0()
    call neighbour_2(2)
    call neighbour_2P(2*p+2)
    call fnf(p)
    deallocate(neighbour_matrix_2)
    deallocate(neighbour_matrix_2P)
    close(10)
    close(20)
    close(30)
    close(40)
end program main

复现

figure 1

module HR
    implicit none
    real,parameter :: h=0.01
    integer,parameter :: N=100,MaxT=360000
    integer,allocatable :: neighbour_matrix_2(:,:),neighbour_matrix_2P(:,:)
    real :: x(N),y(N),z(N),E(N)
contains

    subroutine x0_y0_z0()
        implicit none
        integer :: k
        real :: x1,x2,x3
        call random_seed()
        do k=1,N,1
            call random_number(x1)
            call random_number(x2)
            call random_number(x3)
            x(k)=0.001*(k-N/2.0)+x1
            y(k)=0.002*(k-N/2.0)+x2
            z(k)=0.003*(k-N/2.0)+x3
            write(10,*) x(k),y(k),z(k)
        end do
    end subroutine x0_y0_z0

    subroutine fnf(p)
        implicit none
        integer :: i,j,p,t
        real :: func_x(N),func_y(N),func_z(N),func_e(N),a,b,d,k1,k2,k3,k4,II,JJ(N),C(N),r,s,x_i0,E_ext,x_s,&
                &lambda,thet_s,C1(N),C2(N)
        k1=0.0
        k2=0.0
        k3=0.0
        k4=9.0
        a=1.0
        b=3.0
        d=5.0
        r=0.01
        s=5.0
        x_s=2.0
        x_i0=-1.6
        lambda=10.0
        thet_s=-0.25
        E_ext=0.0
        E=0.0
        II=3.5
        do t=1,MaxT,1
            if(mod(t,5000)==0) then
                write(*,*) t
            end if
            JJ=0.0
            C1=0.0
            C2=0.0
            C=0.0
            do i=1,N,1
                do j=1,N,1
                    JJ(i)=JJ(i)+neighbour_matrix_2(i,j)*(x(j)-x(i))
                    C1(i)=C1(i)+neighbour_matrix_2P(i,j)*(1.0/(1.0+exp(-lambda*(x(j)-thet_s))))
                    C2(i)=C2(i)+neighbour_matrix_2(i,j)*(1.0/(1.0+exp(-lambda*(x(j)-thet_s))))
                end do
                C(i)=C1(i)-C2(i)
                func_x(i)=x(i)+h*(y(i)-a*x(i)**3+b*x(i)**2-z(i)+II+k3*JJ(i)+((-k4/(2.0*p-2.0))*(x_s-x(i))*C(i)))
                func_y(i)=y(i)+h*(1.0-d*x(i)**2-y(i))
                func_z(i)=z(i)+h*(r*(s*(x(i)-x_i0)-z(i)))
                if(mod(t,50)==0.and.t>300000) then
                    write(20,*) i,t-300000,x(i)
                end if
                if(i==5.and.t>300000) then
                    write(30,*) t-300000,x(i)
                    write(70,*) x(i),y(i),z(i)
                    write(*,*) t,x(i)
                end if
                if(i==50.and.t>300000) then
                    write(40,*) t-300000,x(i)
                end if
            end do
            x=func_x
            y=func_y
            z=func_z
        end do
        return
    end subroutine fnf


    !建立网状结构

    !k=2
    subroutine neighbour_2(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix_2(N,N))
        !初始化网络矩阵
        neighbour_matrix_2=0
        do i=1,N-1,1
            do j=i+1,N,1
                if(abs(i-j)<=K/2.or.abs(i-j)>=N-K/2) then
                    neighbour_matrix_2(i,j)=1
                    neighbour_matrix_2(j,i)=1
                end if
            end do
        end do
        !输出矩阵
        do i=1,N,1
            do j=1,N,1
                write(50,"(I2)",advance='no') neighbour_matrix_2(i,j)
            end do
            write(50,*)
        end do
        return
    end subroutine neighbour_2
    !k=2*p+2
    subroutine neighbour_2P(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix_2P(N,N))
        !初始化网络矩阵
        neighbour_matrix_2P=0
        do i=1,N-1,1
            do j=i+1,N,1
                if((abs(i-j)<=K/2.or.abs(i-j)>=N-K/2).and.abs(i-j)/=1) then
                    neighbour_matrix_2P(i,j)=1
                    neighbour_matrix_2P(j,i)=1
                end if
            end do
        end do
        !输出矩阵
        do i=1,N,1
            do j=1,N,1
                write(60,"(I2)",advance='no') neighbour_matrix_2P(i,j)
            end do
            write(60,*)
        end do
        return
    end subroutine neighbour_2P

end module HR



program main
    use HR
    implicit none
    integer :: p
    p=40
    open(10,file="data_x0_y0_z0.txt")
    open(20,file="data_i_t_x.txt")
    open(30,file="data_t_x5.txt")
    open(40,file="data_t_x50.txt")
    open(50,file="neighbour_matrix_2.txt")
    open(60,file="neighbour_matrix_2P.txt")
    open(70,file="x.txt")
    call x0_y0_z0()
    call neighbour_2(2)
    call neighbour_2P(2*p+2)
    call fnf(p)
    deallocate(neighbour_matrix_2)
    deallocate(neighbour_matrix_2P)
    close(10)
    close(20)
    close(30)
    close(40)
    close(50)
    close(60)
    close(70)
end program main

figure 1_A

figure 1_B

figure 2

module HR
    implicit none
    real,parameter :: h=0.01
    integer,parameter :: N=100,MaxT=360000
    integer,allocatable :: neighbour_matrix_2(:,:),neighbour_matrix_2P(:,:)
    real :: x(N),y(N),z(N),E(N)
contains

    subroutine x0_y0_z0()
        implicit none
        integer :: k
        real :: x1,x2,x3
        call random_seed()
        do k=1,N,1
            call random_number(x1)
            call random_number(x2)
            call random_number(x3)
            x(k)=0.001*(k-N/2.0)+x1
            y(k)=0.002*(k-N/2.0)+x2
            z(k)=0.003*(k-N/2.0)+x3
            write(101,*) x(k),y(k),z(k)
        end do
    end subroutine x0_y0_z0

    subroutine fnf(p)
        implicit none
        integer :: i,j,p,t
        real :: func_x(N),func_y(N),func_z(N),func_e(N),a,b,d,k1,k2,k3,k4,II,JJ(N),C(N),r,s,x_i0,E_ext,x_s,&
                &lambda,thet_s,C1(N),C2(N)
        k1=0.0
        k2=0.0
        k3=0.0
        k4=9.0
        a=1.0
        b=3.0
        d=5.0
        r=0.01
        s=5.0
        x_s=2.0
        x_i0=-1.6
        lambda=10.0
        thet_s=-0.25
        E_ext=0.0
        E=0.0
        II=35
        do t=1,MaxT,1
            if(mod(t,5000)==0) then
                write(*,*) t
            end if
            JJ=0.0
            C1=0.0
            C2=0.0
            C=0.0
            do i=1,N,1
                do j=1,N,1
                    JJ(i)=JJ(i)+neighbour_matrix_2(i,j)*(x(j)-x(i))
                    C1(i)=C1(i)+neighbour_matrix_2P(i,j)*(1.0/(1.0+exp(-lambda*(x(j)-thet_s))))
                    C2(i)=C2(i)+neighbour_matrix_2(i,j)*(1.0/(1.0+exp(-lambda*(x(j)-thet_s))))
                end do
                C(i)=C1(i)-C2(i)
                func_x(i)=x(i)+h*(y(i)-a*x(i)**3+b*x(i)**2-z(i)+II+k3*JJ(i)+((-k4/(2.0*p-2.0))*(x_s-x(i))*C(i)))
                func_y(i)=y(i)+h*(1.0-d*x(i)**2-y(i))
                func_z(i)=z(i)+h*(r*(s*(x(i)-x_i0)-z(i)))
                if(t>300000) then
                    write(102,*) i,t-300000,x(i)
                end if
                if(mod(i,2)==0.and.t>300000) then
                    write(i,*) t-300000,x(i)+6.0*i/2.0
                end if
            end do
            x=func_x
            y=func_y
            z=func_z
        end do
        return
    end subroutine fnf


    !建立网状结构

    !k=2
    subroutine neighbour_2(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix_2(N,N))
        !初始化网络矩阵
        neighbour_matrix_2=0
        do i=1,N-1,1
            do j=i+1,N,1
                if(abs(i-j)<=K/2.or.abs(i-j)>=N-K/2) then
                    neighbour_matrix_2(i,j)=1
                    neighbour_matrix_2(j,i)=1
                end if
            end do
        end do
        !输出矩阵
        do i=1,N,1
            do j=1,N,1
                write(103,"(I2)",advance='no') neighbour_matrix_2(i,j)
            end do
            write(103,*)
        end do
        return
    end subroutine neighbour_2
    !k=2*p+2
    subroutine neighbour_2P(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix_2P(N,N))
        !初始化网络矩阵
        neighbour_matrix_2P=0
        do i=1,N-1,1
            do j=i+1,N,1
                if((abs(i-j)<=K/2.or.abs(i-j)>=N-K/2).and.abs(i-j)/=1) then
                    neighbour_matrix_2P(i,j)=1
                    neighbour_matrix_2P(j,i)=1
                end if
            end do
        end do
        !输出矩阵
        do i=1,N,1
            do j=1,N,1
                write(104,"(I2)",advance='no') neighbour_matrix_2P(i,j)
            end do
            write(104,*)
        end do
        return
    end subroutine neighbour_2P

end module HR



program main
    use HR
    implicit none
    integer :: p
    p=40
    open(101,file="data_x0_y0_z0.txt")
    open(102,file="data_i_t_x.txt")
    open(103,file="neighbour_matrix_2.txt")
    open(104,file="neighbour_matrix_2P.txt")
    call x0_y0_z0()
    call neighbour_2(2)
    call neighbour_2P(2*p+2)
    call fnf(p)
    deallocate(neighbour_matrix_2)
    deallocate(neighbour_matrix_2P)
    close(101)
    close(102)
    close(103)
    close(104)
end program main

figure 2_A

figure 2_B

figure 2_C

figure 2_D

figure 3

module HR
    implicit none
    real,parameter :: h=0.01
    integer,parameter :: N=100,MaxT=6079600
    integer,allocatable :: neighbour_matrix_2(:,:),neighbour_matrix_2P(:,:)
    real :: x(N),y(N),z(N),E(N)
contains

    subroutine x0_y0_z0()
        implicit none
        integer :: k
        real :: x1,x2,x3
        call random_seed()
        do k=1,N,1
            !call random_number(x1)
            !call random_number(x2)
            !call random_number(x3)
            !x(k)=0.001*(k-N/2.0)+x1
            !y(k)=0.002*(k-N/2.0)+x2
            !z(k)=0.003*(k-N/2.0)+x3
            read(101,*) x(k),y(k),z(k)
        end do
    end subroutine x0_y0_z0

    subroutine fnf(p)
        implicit none
        integer :: i,j,p,t
        real :: func_x(N),func_y(N),func_z(N),func_e(N),a,b,d,k1,k2,k3,k4,II,JJ(N),C(N),r,s,x_i0,E_ext,x_s,&
                &lambda,thet_s,C1(N),C2(N)
        k1=0.0
        k2=0.0
        k3=0.0
        k4=10.0
        a=1.0
        b=3.0
        d=5.0
        r=0.01
        s=5.0
        x_s=2.0
        x_i0=-1.6
        lambda=10.0
        thet_s=-0.25
        E_ext=0.0
        E=0.0
        II=35
        do t=1,MaxT,1
            if(mod(t,5000)==0) then
                write(*,*) t
            end if
            JJ=0.0
            C1=0.0
            C2=0.0
            C=0.0
            do i=1,N,1
                do j=1,N,1
                    JJ(i)=JJ(i)+neighbour_matrix_2(i,j)*(x(j)-x(i))
                    C1(i)=C1(i)+neighbour_matrix_2P(i,j)*(1.0/(1.0+exp(-lambda*(x(j)-thet_s))))
                    C2(i)=C2(i)+neighbour_matrix_2(i,j)*(1.0/(1.0+exp(-lambda*(x(j)-thet_s))))
                end do
                C(i)=C1(i)-C2(i)
                func_x(i)=x(i)+h*(y(i)-a*x(i)**3+b*x(i)**2-z(i)+II+k3*JJ(i)+((-k4/(2.0*p-2.0))*(x_s-x(i))*C(i)))
                func_y(i)=y(i)+h*(1.0-d*x(i)**2-y(i))
                func_z(i)=z(i)+h*(r*(s*(x(i)-x_i0)-z(i)))
                if(t>5799600.and.mod(t,500)==0) then
                    write(102,*) i,t-5799600,x(i)
                end if
                if(i==60.and.t>6059600) then
                    write(105,*) t-6059600,x(i)
                end if
                if(mod(i,2)==0.and.t>6059600) then
                    write(i,*) t-6059600,x(i)+6.0*(i/2.0)
                end if
            end do
            x=func_x
            y=func_y
            z=func_z
        end do
        return
    end subroutine fnf


    !建立网状结构

    !k=2
    subroutine neighbour_2(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix_2(N,N))
        !初始化网络矩阵
        neighbour_matrix_2=0
        do i=1,N-1,1
            do j=i+1,N,1
                if(abs(i-j)<=K/2.or.abs(i-j)>=N-K/2) then
                    neighbour_matrix_2(i,j)=1
                    neighbour_matrix_2(j,i)=1
                end if
            end do
        end do
        !输出矩阵
        do i=1,N,1
            do j=1,N,1
                write(103,"(I2)",advance='no') neighbour_matrix_2(i,j)
            end do
            write(103,*)
        end do
        return
    end subroutine neighbour_2
    !k=2*p+2
    subroutine neighbour_2P(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix_2P(N,N))
        !初始化网络矩阵
        neighbour_matrix_2P=0
        do i=1,N-1,1
            do j=i+1,N,1
                if((abs(i-j)<=K/2.or.abs(i-j)>=N-K/2).and.abs(i-j)/=1) then
                    neighbour_matrix_2P(i,j)=1
                    neighbour_matrix_2P(j,i)=1
                end if
            end do
        end do
        !输出矩阵
        do i=1,N,1
            do j=1,N,1
                write(104,"(I2)",advance='no') neighbour_matrix_2P(i,j)
            end do
            write(104,*)
        end do
        return
    end subroutine neighbour_2P

end module HR



program main
    use HR
    implicit none
    integer :: p
    p=40
    open(101,file="data_x0_y0_z0.txt")
    open(102,file="data_i_t_x.txt")
    open(103,file="neighbour_matrix_2.txt")
    open(104,file="neighbour_matrix_2P.txt")
    open(105,file="t_x_60.txt")
    call x0_y0_z0()
    call neighbour_2(2)
    call neighbour_2P(2*p+2)
    call fnf(p)
    deallocate(neighbour_matrix_2)
    deallocate(neighbour_matrix_2P)
    close(101)
    close(102)
    close(103)
    close(104)
end program main

figure 3_A

figure 3_B

figure 3_C

figure 3_D

figure 4

module HR
    implicit none
    real,parameter :: h=0.01
    integer,parameter :: N=100,MaxT=90000
    integer,allocatable :: neighbour_matrix_2(:,:),neighbour_matrix_3(:,:),neighbour_matrix_2P(:,:)
    real :: x(N),y(N),z(N),E(N)
contains

    subroutine x0_y0_z0()
        implicit none
        integer :: k
        real :: x1,x2,x3
        call random_seed()
        do k=1,N,1
            call random_number(x1)
            call random_number(x2)
            call random_number(x3)
            x(k)=0.001*(k-N/2.0)+x1
            y(k)=0.002*(k-N/2.0)+x2
            z(k)=0.003*(k-N/2.0)+x3
            write(101,*) x(k),y(k),z(k)
        end do
    end subroutine x0_y0_z0

    subroutine fnf(p)
        implicit none
        integer :: i,j,p,t
        real :: func_x(N),func_y(N),func_z(N),func_e(N),a,b,d,k1,k2,k3,k4,II,JJ(N),C(N),r,s,x_i0,E_ext,x_s,&
                &lambda,thet_s,C1(N),C2(N),Li(N)
        k1=0.0
        k2=0.0
        k3=0.0
        k4=10.0
        a=1.0
        b=3.0
        d=5.0
        r=0.01
        s=5.0
        x_s=2.0
        x_i0=-1.6
        lambda=10.0
        thet_s=-0.25
        E_ext=0.0
        E=0.0
        II=35.0
        do t=1,MaxT,1
            if(mod(t,5000)==0) then
                write(*,*) t
            end if
            JJ=0.0
            C1=0.0
            C2=0.0
            C=0.0
            do i=1,N,1
                do j=1,N,1
                    JJ(i)=JJ(i)+neighbour_matrix_2(i,j)*(x(j)-x(i))
                    C1(i)=C1(i)+neighbour_matrix_2P(i,j)*(1.0/(1.0+exp(-lambda*(x(j)-thet_s))))
                    C2(i)=C2(i)+neighbour_matrix_2(i,j)*(1.0/(1.0+exp(-lambda*(x(j)-thet_s))))
                end do
                C(i)=C1(i)-C2(i)
                func_x(i)=x(i)+h*(y(i)-a*x(i)**3+b*x(i)**2-z(i)+II+k3*JJ(i)+((-k4/(2.0*p-2.0))*(x_s-x(i))*C(i)))
                func_y(i)=y(i)+h*(1.0-d*x(i)**2-y(i))
                func_z(i)=z(i)+h*(r*(s*(x(i)-x_i0)-z(i)))
                if(t>30000.and.mod(t,50)==0) then
                    call parameter_Li(x,y,i,Li)
                    write(102,*) i,t-30000,Li(i)
                end if
            end do
            x=func_x
            y=func_y
            z=func_z
        end do
        return
    end subroutine fnf

    subroutine parameter_Li(x,y,i,Li)
        implicit none
        real :: x(N),y(N),Ans_L_X,Ans_L_Y,Li(N),eta
        integer :: i,j
        Ans_L_X=0.0
        Ans_L_Y=0.0
        eta=2.0
        do j=1,N,1
            Ans_L_X=Ans_L_X+neighbour_matrix_3(i,j)*(cos(atan(y(j)/x(j))))
            Ans_L_Y=Ans_L_Y+neighbour_matrix_3(i,j)*(sin(atan(y(j)/x(j))))
        end do
        Li(i)=sqrt((Ans_L_X/(2.0*eta))**2+(Ans_L_Y/(2.0*eta))**2)
        return
    end subroutine parameter_Li

    !建立网状结构

    !k=2
    subroutine neighbour_2(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix_2(N,N))
        !初始化网络矩阵
        neighbour_matrix_2=0
        do i=1,N-1,1
            do j=i+1,N,1
                if(abs(i-j)<=K/2.or.abs(i-j)>=N-K/2) then
                    neighbour_matrix_2(i,j)=1
                    neighbour_matrix_2(j,i)=1
                end if
            end do
        end do
        !输出矩阵
        do i=1,N,1
            do j=1,N,1
                write(103,"(I2)",advance='no') neighbour_matrix_2(i,j)
            end do
            write(103,*)
        end do
        return
    end subroutine neighbour_2
    !k=2
    subroutine neighbour_3(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix_3(N,N))
        !初始化网络矩阵
        neighbour_matrix_3=0
        do i=1,N-1,1
            do j=i+1,N,1
                if((abs(i-j)<=K/2.or.abs(i-j)>=N-K/2).and.abs(i-j)<=2) then
                    neighbour_matrix_3(i,j)=1
                    neighbour_matrix_3(j,i)=1
                end if
            end do
        end do
        !输出矩阵
        do i=1,N,1
            do j=1,N,1
                write(105,"(I2)",advance='no') neighbour_matrix_3(i,j)
            end do
            write(105,*)
        end do
        return
    end subroutine neighbour_3
    !k=2*p+2
    subroutine neighbour_2P(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix_2P(N,N))
        !初始化网络矩阵
        neighbour_matrix_2P=0
        do i=1,N-1,1
            do j=i+1,N,1
                if((abs(i-j)<=K/2.or.abs(i-j)>=N-K/2).and.abs(i-j)/=1) then
                    neighbour_matrix_2P(i,j)=1
                    neighbour_matrix_2P(j,i)=1
                end if
            end do
        end do
        !输出矩阵
        do i=1,N,1
            do j=1,N,1
                write(104,"(I2)",advance='no') neighbour_matrix_2P(i,j)
            end do
            write(104,*)
        end do
        return
    end subroutine neighbour_2P

end module HR



program main
    use HR
    implicit none
    integer :: p
    p=40
    open(101,file="data_x0_y0_z0.txt")
    open(102,file="data_i_t_Li.txt")
    open(103,file="neighbour_matrix_2.txt")
    open(104,file="neighbour_matrix_2P.txt")
    open(105,file="neighbour_matrix_3.txt")
    call x0_y0_z0()
    call neighbour_2(2)
    call neighbour_3(4)
    call neighbour_2P(2*p+2)
    call fnf(p)
    deallocate(neighbour_matrix_2)
    deallocate(neighbour_matrix_2P)
    close(101)
    close(102)
    close(103)
    close(104)
    close(105)
end program main

figure 4_A

figure 4_B

figure 5

module HR
    implicit none
    real,parameter :: h=0.01,PI=3.1415926
    integer,parameter :: N=100,MaxT=260000
    integer,allocatable :: neighbour_matrix_2(:,:),neighbour_matrix_2P(:,:)
    real :: x(N),y(N),z(N),E(N)
contains

    subroutine x0_y0_z0()
        implicit none
        integer :: k
        real :: x1,x2,x3
        call random_seed()
        do k=1,N,1
            call random_number(x1)
            call random_number(x2)
            call random_number(x3)
            x(k)=0.001*(k-N/2.0)+x1
            y(k)=0.002*(k-N/2.0)+x2
            z(k)=0.003*(k-N/2.0)+x3
            write(101,*) x(k),y(k),z(k)
        end do
    end subroutine x0_y0_z0

    subroutine fnf(p)
        implicit none
        integer :: i,j,p,t
        real :: func_x(N),func_y(N),func_z(N),func_e(N),a,b,d,k1,k2,k3,k4,II,JJ(N),C(N),r,s,x_i0,E_ext,x_s,&
                &lambda,thet_s,C1(N),C2(N),Li(N),Em,f
        k1=0.7
        k2=0.001
        k3=0.0
        k4=9.0
        a=1.0
        b=3.0
        d=5.0
        r=0.01
        s=5.0
        x_s=2.0
        x_i0=-1.6
        lambda=10.0
        thet_s=-0.25
        E_ext=0.0
        E=0.0
        Em=1.5
        f=0.01
        II=3.5
        do t=1,MaxT,1
            if(mod(t,5000)==0) then
                write(*,*) t
            end if
            JJ=0.0
            C1=0.0
            C2=0.0
            C=0.0
            do i=1,N,1
                do j=1,N,1
                    JJ(i)=JJ(i)+neighbour_matrix_2(i,j)*(x(j)-x(i))
                    C1(i)=C1(i)+neighbour_matrix_2P(i,j)*(1.0/(1.0+exp(-lambda*(x(j)-thet_s))))
                    C2(i)=C2(i)+neighbour_matrix_2(i,j)*(1.0/(1.0+exp(-lambda*(x(j)-thet_s))))
                end do
                C(i)=C1(i)-C2(i)
                func_x(i)=x(i)+h*(y(i)-a*x(i)**3+b*x(i)**2-z(i)+II+k3*JJ(i)+((-k4/(2.0*p-2.0))*(x_s-x(i))*C(i)))
                func_y(i)=y(i)+h*(1.0-d*x(i)**2-y(i)+k1*E(i))
                func_z(i)=z(i)+h*(r*(s*(x(i)-x_i0)-z(i)))
                if(i>50.and.t>100000) then
                    func_e(i)=E(i)+h*(k2*y(i)+Em*sin(2.0*PI*f*(t-100000)*h))
                else
                    func_e(i)=0.0
                end if
                if(t>200000.and.mod(t,50)==0) then
                    write(102,*) i,t-200000,x(i)
                end if
                if(t>200000.and.i==5) then
                    write(105,*) t-200000,x(i)
                end if   
                if(t>200000.and.i==55) then
                    write(106,*) t-200000,x(i)
                end if
            end do
            x=func_x
            y=func_y
            z=func_z
            E=func_e
        end do
        return
    end subroutine fnf

    !建立网状结构

    !k=2
    subroutine neighbour_2(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix_2(N,N))
        !初始化网络矩阵
        neighbour_matrix_2=0
        do i=1,N-1,1
            do j=i+1,N,1
                if(abs(i-j)<=K/2.or.abs(i-j)>=N-K/2) then
                    neighbour_matrix_2(i,j)=1
                    neighbour_matrix_2(j,i)=1
                end if
            end do
        end do
        !输出矩阵
        do i=1,N,1
            do j=1,N,1
                write(103,"(I2)",advance='no') neighbour_matrix_2(i,j)
            end do
            write(103,*)
        end do
        return
    end subroutine neighbour_2

    !k=2*p+2
    subroutine neighbour_2P(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix_2P(N,N))
        !初始化网络矩阵
        neighbour_matrix_2P=0
        do i=1,N-1,1
            do j=i+1,N,1
                if((abs(i-j)<=K/2.or.abs(i-j)>=N-K/2).and.abs(i-j)/=1) then
                    neighbour_matrix_2P(i,j)=1
                    neighbour_matrix_2P(j,i)=1
                end if
            end do
        end do
        !输出矩阵
        do i=1,N,1
            do j=1,N,1
                write(104,"(I2)",advance='no') neighbour_matrix_2P(i,j)
            end do
            write(104,*)
        end do
        return
    end subroutine neighbour_2P

end module HR



program main
    use HR
    implicit none
    integer :: p
    p=40
    open(101,file="data_x0_y0_z0.txt")
    open(102,file="data_i_t_x.txt")
    open(103,file="neighbour_matrix_2.txt")
    open(104,file="neighbour_matrix_2P.txt")
    open(105,file="t_x_5.txt")
    open(106,file="t_x_55.txt")
    call x0_y0_z0()
    call neighbour_2(2)
    call neighbour_2P(2*p+2)
    call fnf(p)
    deallocate(neighbour_matrix_2)
    deallocate(neighbour_matrix_2P)
    close(101)
    close(102)
    close(103)
    close(104)
end program main

figure 5_A

figure 5_B

figure 5_C

module HR
    implicit none
    real,parameter :: h=0.01,PI=3.1415926
    integer,parameter :: N=100,MaxT=160000
    integer,allocatable :: neighbour_matrix_2(:,:),neighbour_matrix_2P(:,:)
    real :: x(N),y(N),z(N),E(N)
contains

    subroutine x0_y0_z0()
        implicit none
        integer :: k
        real :: x1,x2,x3
        call random_seed()
        do k=1,N,1
            call random_number(x1)
            call random_number(x2)
            call random_number(x3)
            x(k)=0.001*(k-N/2.0)+x1
            y(k)=0.002*(k-N/2.0)+x2
            z(k)=0.003*(k-N/2.0)+x3
            write(101,*) x(k),y(k),z(k)
        end do
    end subroutine x0_y0_z0

    subroutine fnf(p)
        implicit none
        integer :: i,j,p,t
        real :: func_x(N),func_y(N),func_z(N),func_e(N),a,b,d,k1,k2,k3,k4,II,JJ(N),C(N),r,s,x_i0,E_ext,x_s,&
                &lambda,thet_s,C1(N),C2(N),Li(N),Em,f
        k1=0.7
        k2=0.001
        k3=0.0
        k4=9.0
        a=1.0
        b=3.0
        d=5.0
        r=0.01
        s=5.0
        x_s=2.0
        x_i0=-1.6
        lambda=10.0
        thet_s=-0.25
        E_ext=0.0
        E=0.0
        Em=1.5
        f=12.0
        II=3.5
        do t=1,MaxT,1
            if(mod(t,5000)==0) then
                write(*,*) t
            end if
            JJ=0.0
            C1=0.0
            C2=0.0
            C=0.0
            do i=1,N,1
                do j=1,N,1
                    JJ(i)=JJ(i)+neighbour_matrix_2(i,j)*(x(j)-x(i))
                    C1(i)=C1(i)+neighbour_matrix_2P(i,j)*(1.0/(1.0+exp(-lambda*(x(j)-thet_s))))
                    C2(i)=C2(i)+neighbour_matrix_2(i,j)*(1.0/(1.0+exp(-lambda*(x(j)-thet_s))))
                end do
                C(i)=C1(i)-C2(i)
                func_x(i)=x(i)+h*(y(i)-a*x(i)**3+b*x(i)**2-z(i)+II+k3*JJ(i)+((-k4/(2.0*p-2.0))*(x_s-x(i))*C(i)))
                func_y(i)=y(i)+h*(1.0-d*x(i)**2-y(i)+k1*E(i))
                func_z(i)=z(i)+h*(r*(s*(x(i)-x_i0)-z(i)))
                if(i>50) then
                    func_e(i)=E(i)+h*(k2*y(i)+Em*sin(2.0*PI*f*t*h))
                else
                    func_e(i)=0.0
                end if
                if(t>100000.and.mod(t,200)==0) then
                    write(102,*) i,t-100000,x(i)
                end if
                if(t>100000.and.i==5) then
                    write(105,*) t-100000,x(i)
                end if   
                if(t>100000.and.i==55) then
                    write(106,*) t-100000,x(i)
                end if
            end do
            x=func_x
            y=func_y
            z=func_z
            E=func_e
        end do
        return
    end subroutine fnf

    !建立网状结构

    !k=2
    subroutine neighbour_2(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix_2(N,N))
        !初始化网络矩阵
        neighbour_matrix_2=0
        do i=1,N-1,1
            do j=i+1,N,1
                if(abs(i-j)<=K/2.or.abs(i-j)>=N-K/2) then
                    neighbour_matrix_2(i,j)=1
                    neighbour_matrix_2(j,i)=1
                end if
            end do
        end do
        !输出矩阵
        do i=1,N,1
            do j=1,N,1
                write(103,"(I2)",advance='no') neighbour_matrix_2(i,j)
            end do
            write(103,*)
        end do
        return
    end subroutine neighbour_2

    !k=2*p+2
    subroutine neighbour_2P(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix_2P(N,N))
        !初始化网络矩阵
        neighbour_matrix_2P=0
        do i=1,N-1,1
            do j=i+1,N,1
                if((abs(i-j)<=K/2.or.abs(i-j)>=N-K/2).and.abs(i-j)/=1) then
                    neighbour_matrix_2P(i,j)=1
                    neighbour_matrix_2P(j,i)=1
                end if
            end do
        end do
        !输出矩阵
        do i=1,N,1
            do j=1,N,1
                write(104,"(I2)",advance='no') neighbour_matrix_2P(i,j)
            end do
            write(104,*)
        end do
        return
    end subroutine neighbour_2P

end module HR



program main
    use HR
    implicit none
    integer :: p
    p=40
    open(101,file="data_x0_y0_z0.txt")
    open(102,file="data_i_t_x.txt")
    open(103,file="neighbour_matrix_2.txt")
    open(104,file="neighbour_matrix_2P.txt")
    open(105,file="t_x_5.txt")
    open(106,file="t_x_55.txt")
    call x0_y0_z0()
    call neighbour_2(2)
    call neighbour_2P(2*p+2)
    call fnf(p)
    deallocate(neighbour_matrix_2)
    deallocate(neighbour_matrix_2P)
    close(101)
    close(102)
    close(103)
    close(104)
end program main

figure 5_D

figure 6

figure 6_A

module HR
    implicit none
    real,parameter :: h=0.01,PI=3.1415926
    integer,parameter :: N=100,MaxT=160000
    integer,allocatable :: neighbour_matrix_2(:,:),neighbour_matrix_2P(:,:)
    real :: x(N),y(N),z(N),E(N)
contains

    subroutine x0_y0_z0()
        implicit none
        integer :: k
        real :: x1,x2,x3
        call random_seed()
        do k=1,N,1
            call random_number(x1)
            call random_number(x2)
            call random_number(x3)
            x(k)=0.001*(k-N/2.0)+x1
            y(k)=0.002*(k-N/2.0)+x2
            z(k)=0.003*(k-N/2.0)+x3
            write(101,*) x(k),y(k),z(k)
        end do
    end subroutine x0_y0_z0

    subroutine fnf(p)
        implicit none
        integer :: i,j,p,t
        real :: func_x(N),func_y(N),func_z(N),func_e(N),a,b,d,k1,k2,k3,k4,II,JJ(N),C(N),r,s,x_i0,E_ext,x_s,&
                &lambda,thet_s,C1(N),C2(N),Li(N),Em,f
        k1=0.7
        k2=0.001
        k3=0.0
        k4=9.0
        a=1.0
        b=3.0
        d=5.0
        r=0.01
        s=5.0
        x_s=2.0
        x_i0=-1.6
        lambda=10.0
        thet_s=-0.25
        E_ext=0.0
        E=0.0
        Em=1.5
        f=12.0
        II=3.5
        do t=1,MaxT,1
            if(mod(t,5000)==0) then
                write(*,*) t
            end if
            JJ=0.0
            C1=0.0
            C2=0.0
            C=0.0
            do i=1,N,1
                do j=1,N,1
                    JJ(i)=JJ(i)+neighbour_matrix_2(i,j)*(x(j)-x(i))
                    C1(i)=C1(i)+neighbour_matrix_2P(i,j)*(1.0/(1.0+exp(-lambda*(x(j)-thet_s))))
                    C2(i)=C2(i)+neighbour_matrix_2(i,j)*(1.0/(1.0+exp(-lambda*(x(j)-thet_s))))
                end do
                C(i)=C1(i)-C2(i)
                func_x(i)=x(i)+h*(y(i)-a*x(i)**3+b*x(i)**2-z(i)+II+k3*JJ(i)+((-k4/(2.0*p-2.0))*(x_s-x(i))*C(i)))
                func_y(i)=y(i)+h*(1.0-d*x(i)**2-y(i)+k1*E(i))
                func_z(i)=z(i)+h*(r*(s*(x(i)-x_i0)-z(i)))
                if(i>25) then
                    func_e(i)=E(i)+h*(k2*y(i)+Em*sin(2.0*PI*f*t*h))
                else
                    func_e(i)=0.0
                end if
                if(t>100000.and.mod(t,200)==0) then
                    write(102,*) i,t-100000,x(i)
                end if
            end do
            x=func_x
            y=func_y
            z=func_z
            E=func_e
        end do
        return
    end subroutine fnf

    !建立网状结构

    !k=2
    subroutine neighbour_2(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix_2(N,N))
        !初始化网络矩阵
        neighbour_matrix_2=0
        do i=1,N-1,1
            do j=i+1,N,1
                if(abs(i-j)<=K/2.or.abs(i-j)>=N-K/2) then
                    neighbour_matrix_2(i,j)=1
                    neighbour_matrix_2(j,i)=1
                end if
            end do
        end do
        !输出矩阵
        do i=1,N,1
            do j=1,N,1
                write(103,"(I2)",advance='no') neighbour_matrix_2(i,j)
            end do
            write(103,*)
        end do
        return
    end subroutine neighbour_2

    !k=2*p+2
    subroutine neighbour_2P(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix_2P(N,N))
        !初始化网络矩阵
        neighbour_matrix_2P=0
        do i=1,N-1,1
            do j=i+1,N,1
                if((abs(i-j)<=K/2.or.abs(i-j)>=N-K/2).and.abs(i-j)/=1) then
                    neighbour_matrix_2P(i,j)=1
                    neighbour_matrix_2P(j,i)=1
                end if
            end do
        end do
        !输出矩阵
        do i=1,N,1
            do j=1,N,1
                write(104,"(I2)",advance='no') neighbour_matrix_2P(i,j)
            end do
            write(104,*)
        end do
        return
    end subroutine neighbour_2P

end module HR



program main
    use HR
    implicit none
    integer :: p
    p=40
    open(101,file="data_x0_y0_z0.txt")
    open(102,file="data_i_t_x.txt")
    open(103,file="neighbour_matrix_2.txt")
    open(104,file="neighbour_matrix_2P.txt")
    call x0_y0_z0()
    call neighbour_2(2)
    call neighbour_2P(2*p+2)
    call fnf(p)
    deallocate(neighbour_matrix_2)
    deallocate(neighbour_matrix_2P)
    close(101)
    close(102)
    close(103)
    close(104)
end program main

figure 6_B

figure 6_C

module HR
    implicit none
    real,parameter :: h=0.01,PI=3.1415926
    integer,parameter :: N=100,MaxT=160000
    integer,allocatable :: neighbour_matrix_2(:,:),neighbour_matrix_2P(:,:)
    real(kind=8) :: x(N),y(N),z(N),E(N)
    integer :: Nb
contains

    subroutine x0_y0_z0()
        implicit none
        integer :: k
        real :: x1,x2,x3
        call random_seed()
        do k=1,N,1
            call random_number(x1)
            call random_number(x2)
            call random_number(x3)
            x(k)=0.001*(k-N/2.0)+x1
            y(k)=0.002*(k-N/2.0)+x2
            z(k)=0.003*(k-N/2.0)+x3
!            write(101,*) x(k),y(k),z(k)
        end do
    end subroutine x0_y0_z0

    subroutine fnf(p)
        implicit none
        integer :: i,j,p,t
        real(kind=8) :: func_x(N),func_y(N),func_z(N),func_e(N),a,b,d,k1,k2,k3,k4,II,JJ(N),C(N),r,s,x_i0,E_ext,x_s,&
                &lambda,thet_s,C1(N),C2(N),Li(N),Em,f,xx(N,60000)
        k1=0.7
        k2=0.001
        k3=0.0
        k4=9.0
        a=1.0
        b=3.0
        d=5.0
        r=0.01
        s=5.0
        x_s=2.0
        x_i0=-1.6
        lambda=10.0
        thet_s=-0.25
        E_ext=0.0
        E=0.0
        Em=1.5
        f=12.0
        II=3.5
        do Nb=1,N,1
            xx=0.0
            call x0_y0_z0()
            do t=1,MaxT,1
                JJ=0.0
                C1=0.0
                C2=0.0
                C=0.0
                do i=1,N,1
                    do j=1,N,1
                        JJ(i)=JJ(i)+neighbour_matrix_2(i,j)*(x(j)-x(i))
                        C1(i)=C1(i)+neighbour_matrix_2P(i,j)*(1.0/(1.0+exp(-lambda*(x(j)-thet_s))))
                        C2(i)=C2(i)+neighbour_matrix_2(i,j)*(1.0/(1.0+exp(-lambda*(x(j)-thet_s))))
                    end do
                    C(i)=C1(i)-C2(i)
                    func_x(i)=x(i)+h*(y(i)-a*x(i)**3+b*x(i)**2-z(i)+II+k3*JJ(i)+((-k4/(2.0*p-2.0))*(x_s-x(i))*C(i)))
                    func_y(i)=y(i)+h*(1.0-d*x(i)**2-y(i)+k1*E(i))
                    func_z(i)=z(i)+h*(r*(s*(x(i)-x_i0)-z(i)))
                    if(i>=N-Nb) then
                        func_e(i)=E(i)+h*(k2*y(i)+Em*sin(2.0*PI*f*t*h))
                    else
                        func_e(i)=0.0
                    end if
                    if(t>100000) then
                        xx(i,t-100000)=x(i)
                    end if
                end do
                x=func_x
                y=func_y
                z=func_z
                E=func_e
            end do
            call SIDM(xx)
        end do
        return
    end subroutine fnf

    !SI不相干强度
    subroutine SIDM(xx)
        implicit none
        integer :: i,j,t,m
        real(kind=8) :: xx(N,60000),delta
        real(kind=8) :: sig(Nb),S(Nb),ZZm(N)
        real(kind=8) :: Sm,SI,DM,ZZ,nn
        delta=0.02*abs(maxval(xx(:,:))-minval(xx(:,:)))
        nn=real(N)/real(Nb)
        SI=0.0
        DM=0.0
        ZZm=0.0
        do i=1,N-1,1
            do t=1,60000,1
                ZZm(i)=ZZm(i)+(xx(i,t)-xx(i+1,t))
            end do
        end do
        do t=1,60000,1
            ZZm(N)=ZZm(N)+(xx(N,t)-xx(1,t))
        end do
        ZZm=ZZm/60000.0
        sig=0.0
        do m=1,Nb,1
            do t=1,60000,1
                ZZ=0.0
                do j=nn*(m-1)+1,m*nn,1
                    if(j/=N) then
                        ZZ=ZZ+((xx(j,t)-xx(j+1,t))-ZZm(m))**2
                    else
                        ZZ=ZZ+((xx(N,t)-xx(1,t))-ZZm(m))**2
                    end if
                end do
                sig(m)=sig(m)+sqrt((1.0/nn)*ZZ)
            end do
        end do
        sig=sig/60000.0
        do m=1,Nb,1
            Sm=delta-sig(m)
            write(*,*) sig(m),delta
            if(Sm>0.0) then
                S(m)=1.0
            else
                S(m)=0.0
            end if
        end do
        SI=1.0-(1.0/Nb)*sum(S(:))
        do m=1,Nb-1,1
            DM=DM+abs(S(m)-S(m+1))
        end do
        if(Nb==N) then
            DM=DM+abs(S(N)-S(1))
        end if
        DM=DM/2.0
        write(105,*) Nb,SI
        write(106,*) Nb,DM
        write(*,*) Nb,SI,DM
        return
    end subroutine SIDM

    !建立网状结构

    !k=2
    subroutine neighbour_2(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix_2(N,N))
        !初始化网络矩阵
        neighbour_matrix_2=0
        do i=1,N-1,1
            do j=i+1,N,1
                if(abs(i-j)<=K/2.or.abs(i-j)>=N-K/2) then
                    neighbour_matrix_2(i,j)=1
                    neighbour_matrix_2(j,i)=1
                end if
            end do
        end do
        !输出矩阵
        do i=1,N,1
            do j=1,N,1
                write(103,"(I2)",advance='no') neighbour_matrix_2(i,j)
            end do
            write(103,*)
        end do
        return
    end subroutine neighbour_2

    !k=2*p+2
    subroutine neighbour_2P(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix_2P(N,N))
        !初始化网络矩阵
        neighbour_matrix_2P=0
        do i=1,N-1,1
            do j=i+1,N,1
                if((abs(i-j)<=K/2.or.abs(i-j)>=N-K/2).and.abs(i-j)/=1) then
                    neighbour_matrix_2P(i,j)=1
                    neighbour_matrix_2P(j,i)=1
                end if
            end do
        end do
        !输出矩阵
        do i=1,N,1
            do j=1,N,1
                write(104,"(I2)",advance='no') neighbour_matrix_2P(i,j)
            end do
            write(104,*)
        end do
        return
    end subroutine neighbour_2P

end module HR



program main
    use HR
    implicit none
    integer :: p
    p=40
    open(101,file="data_x0_y0_z0.txt")
    open(102,file="data_i_t_x.txt")
    open(103,file="neighbour_matrix_2.txt")
    open(104,file="neighbour_matrix_2P.txt")
    open(105,file="SI.txt")
    open(106,file="DM.txt")
    open(107,file="Nb_11.txt")
    call neighbour_2(2)
    call neighbour_2P(2*p+2)
    call fnf(p)
    deallocate(neighbour_matrix_2)
    deallocate(neighbour_matrix_2P)
    close(101)
    close(102)
    close(103)
    close(104)
    close(105)
    close(106)
end program main

figure 6_D

    module mini
    implicit none
    integer(kind=4),parameter:: N=100,M=240000
    real(kind=8),parameter::h=0.01,Pi=3.1415926d0
    integer(kind=4)::Nb
    end module mini

    program main
    use mini
    implicit none
    integer(kind=4)::i,j,t
    integer(kind=4),parameter::P=40
    real(kind=8),parameter::a=1.0,b=3.0,d=5.0,r=0.01,s=5.0,xi0=-1.6,thetas=-0.25,lamda=10.0,xs=2.0
    real(kind=8),dimension(1:N)::x,y,z,E,x1,y1,z1,E1
    real(kind=8):: k3,k4,Ie,Ci,Ci1,Ci2,Ji,r1,r2,r3,f,Em,k1,k2
    real(kind=8)::xx(1:N,1:M)


    call random_seed()

    open(unit=11,file='SI.dat')
    open(unit=12,file='DM.dat')

    do Nb=1,100
 !   open(unit=1,file='时空斑图.dat')

    f=12;Em=1.5
    k1=0.7
    k2=0.001
    k3=0
    k4=9
    Ie=3.5

    x=0;y=0;z=0
    do i=1,N
        call random_number(r1)
        call random_number(r2)
        call random_number(r3)
        x(i)=0.001*(i-N/2.0)+r1
        y(i)=0.002*(i-N/2.0)+r2
        z(i)=0.003*(i-N/2.0)+r3
    end do
    E=0
    do t=1,M
        !if(mod(t,10000)==0) print*,t*h  
        do i=1,N
            !-----电突触耦合---------对应（2）式
            Ji=0;
            do j=1,N
                if(abs(i-j)<=1 .or. abs(i-j)>=N-1)then
                    Ji=Ji+(x(j)-x(i))
                end if
            end do

            !-----化学突触耦合---------对应（3）（4）式
            Ci1=0;Ci2=0;Ci=0
            do j=1,N
                if((abs(i-j)<=P .or. abs(i-j)>=N-P) .and. abs(i-j)>1 )then
                    Ci1=Ci1+(1.0/(1.0+exp(-lamda*(x(j)-thetas))))
                end if
                if((abs(i-j)<=1 .or. abs(i-j)>=N-1) .and. abs(i-j)/=0)then
                    Ci2=Ci2+(1.0/(1.0+exp(-lamda*(x(j)-thetas))))
                end if
            end do
            Ci=Ci1-Ci2

            x1(i)=x(i)+h*(y(i) -a*x(i)**3 + b*x(i)**2 - z(i)+ Ie+ k3*Ji + ((-k4/(2.0*p-2.0)*(xs-x(i)))*Ci))
            y1(i)=y(i)+h*(1.0 - d*x(i)**2 - y(i)+k1*E(i))
            z1(i)=z(i)+h*(r*(s*(x(i)-xi0) - z(i)))
            if(i>=N-Nb)then
                E1(i)=E(i)+h*(k2*y(i)+Em*sin(2*Pi*f*t*h))
            else
                E1(i)=0
            end if
            xx(i,t)=x1(i)
        end do

        x=x1;y=y1;z=z1;E=E1
    end do

    call SIDM(xx)

    end do
    stop
    end

subroutine SIDM(xx)
use mini
implicit none
integer(kind=4)::i,j,t,mm
real(kind=8)::Delta!,parameter::delta=0.15!1
real(kind=8)::xx(1:N,1:M),sig(1:Nb),SS(1:Nb),ZZa(1:N)
real(kind=8)::nn,sigmam,Sm,SI,DM,SDM,ZZJ

!print*,"SI,DM",0.02*abs(maxval(xx(:,:))-minval(xx(:,:)))
delta=0.02*abs(maxval(xx(:,:))-minval(xx(:,:)))
nn=real(N)/real(Nb)
SI=0
!--------<ZZm> 算法--------
ZZa=0
do i=1,N-1
    do t=M/2+1,M
        Zza(i)=Zza(i)+xx(i,t)-xx(i+1,t)
    end do
end do
do t=M/2+1,M
    Zza(N)=Zza(N)+xx(N,t)-xx(1,t)
end do
Zza=zza/real(M/2)  
!--------<ZZm> 算法--------
!-----------------------sigma(m)的算法------------------------
sig=0
do mm=1,Nb
    !给定确定的m值，先算出sigma(m)的值，即m确定时sigma(m,t)关于时间的平均
    do t=M/2+1,M
        zzj=0   !计算 (zzj-<zzm>)^2
        do j=nn*(mm-1)+1,nn*mm
            if(j/=N)then
                zzj=zzj+(xx(j,t)-xx(j+1,t)-zza(mm))**2
            else
                zzj=zzj+(xx(N,t)-xx(1,t)-zza(mm))**2
            end if
        end do
        sig(mm)=sig(mm)+sqrt(1.0/nn*zzj)
    end do
end do
sig=sig/real(M/2)
!-----------------------sigma(m)的算法------------------------

!-----------------------Sm的计算方法-------------------
do mm=1,Nb
    Sm=delta-sig(mm)
    write(*,*) sig(mm),delta
    if(Sm>0)then
        SS(mm)=1
    else
        SS(mm)=0
    end if
end do
!-----------------------Sm的计算方法-------------------
SI=1-1.0/Nb*sum(SS(:))

DM=0
SDM=0
do mm=1,Nb-1
    SDM=SDM+abs(ss(mm)-ss(mm+1))
end do
SDM=SDM+abs(ss(Nb)-ss(1))
DM=SDM/2.0

write(11,*)Nb,SI
write(12,*)Nb,DM
print*,Nb,SI,DM

return
end subroutine SIDM

figure 7

module HR
    implicit none
    real,parameter :: h=0.01,PI=3.1415926
    integer,parameter :: N=100,MaxT=160000
    integer,allocatable :: neighbour_matrix_2(:,:),neighbour_matrix_2P(:,:)
    real :: x(N),y(N),z(N),E(N)
contains

    subroutine x0_y0_z0()
        implicit none
        integer :: k
        real :: x1,x2,x3
        call random_seed()
        do k=1,N,1
            call random_number(x1)
            call random_number(x2)
            call random_number(x3)
            x(k)=0.001*(k-N/2.0)+x1
            y(k)=0.002*(k-N/2.0)+x2
            z(k)=0.003*(k-N/2.0)+x3
            read(101,*) x(k),y(k),z(k)
        end do
    end subroutine x0_y0_z0

    subroutine fnf(p)
        implicit none
        integer :: i,j,p,t
        real :: func_x(N),func_y(N),func_z(N),func_e(N),a,b,d,k1,k2,k3,k4,II,JJ(N),C(N),r,s,x_i0,E_ext,x_s,&
                &lambda,thet_s,C1(N),C2(N),Li(N),Em,f
        k1=0.7
        k2=0.001
        k3=0.0
        k4=9.0
        a=1.0
        b=3.0
        d=5.0
        r=0.01
        s=5.0
        x_s=2.0
        x_i0=-1.6
        lambda=10.0
        thet_s=-0.25
        E_ext=0.0
        E=0.0
        Em=1.5
        f=12.0
        II=3.5
        do t=1,MaxT,1
            if(mod(t,5000)==0) then
                write(*,*) t
            end if
            JJ=0.0
            C1=0.0
            C2=0.0
            C=0.0
            do i=1,N,1
                do j=1,N,1
                    JJ(i)=JJ(i)+neighbour_matrix_2(i,j)*(x(j)-x(i))
                    C1(i)=C1(i)+neighbour_matrix_2P(i,j)*(1.0/(1.0+exp(-lambda*(x(j)-thet_s))))
                    C2(i)=C2(i)+neighbour_matrix_2(i,j)*(1.0/(1.0+exp(-lambda*(x(j)-thet_s))))
                end do
                C(i)=C1(i)-C2(i)
                func_x(i)=x(i)+h*(y(i)-a*x(i)**3+b*x(i)**2-z(i)+II+k3*JJ(i)+((-k4/(2.0*p-2.0))*(x_s-x(i))*C(i)))
                func_y(i)=y(i)+h*(1.0-d*x(i)**2-y(i)+k1*E(i))
                func_z(i)=z(i)+h*(r*(s*(x(i)-x_i0)-z(i)))
                if(i<25.or.(i>50.and.i<75)) then
                    func_e(i)=E(i)+h*(k2*y(i)+Em*sin(2.0*PI*f*t*h))
                else
                    func_e(i)=0.0
                end if
                if(t>100000.and.mod(t,200)==0) then
                    write(102,*) i,t-100000,x(i)
                end if
                if(t==160000) then
                    write(105,*) i,x(i)
                end if
            end do
            x=func_x
            y=func_y
            z=func_z
            E=func_e
        end do
        return
    end subroutine fnf

    !建立网状结构

    !k=2
    subroutine neighbour_2(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix_2(N,N))
        !初始化网络矩阵
        neighbour_matrix_2=0
        do i=1,N-1,1
            do j=i+1,N,1
                if(abs(i-j)<=K/2.or.abs(i-j)>=N-K/2) then
                    neighbour_matrix_2(i,j)=1
                    neighbour_matrix_2(j,i)=1
                end if
            end do
        end do
        !输出矩阵
        do i=1,N,1
            do j=1,N,1
                write(103,"(I2)",advance='no') neighbour_matrix_2(i,j)
            end do
            write(103,*)
        end do
        return
    end subroutine neighbour_2

    !k=2*p+2
    subroutine neighbour_2P(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix_2P(N,N))
        !初始化网络矩阵
        neighbour_matrix_2P=0
        do i=1,N-1,1
            do j=i+1,N,1
                if((abs(i-j)<=K/2.or.abs(i-j)>=N-K/2).and.abs(i-j)/=1) then
                    neighbour_matrix_2P(i,j)=1
                    neighbour_matrix_2P(j,i)=1
                end if
            end do
        end do
        !输出矩阵
        do i=1,N,1
            do j=1,N,1
                write(104,"(I2)",advance='no') neighbour_matrix_2P(i,j)
            end do
            write(104,*)
        end do
        return
    end subroutine neighbour_2P

end module HR



program main
    use HR
    implicit none
    integer :: p
    p=40
    open(101,file="data_x0_y0_z0.txt")
    open(102,file="data_i_t_x.txt")
    open(103,file="neighbour_matrix_2.txt")
    open(104,file="neighbour_matrix_2P.txt")
    open(105,file="i_x.txt")
    call x0_y0_z0()
    call neighbour_2(2)
    call neighbour_2P(2*p+2)
    call fnf(p)
    deallocate(neighbour_matrix_2)
    deallocate(neighbour_matrix_2P)
    close(101)
    close(102)
    close(103)
    close(104)
end program main

figure 7_A

figure 7_B

figure 8

module HR
    implicit none
    real,parameter :: h=0.01,PI=3.1415926
    integer,parameter :: N=20,MaxT=160000
    integer,allocatable :: neighbour_matrix_2(:,:),neighbour_matrix_2P(:,:)
    real(kind=8) :: x(N),y(N),z(N),E(N),xx(N,MaxT)
    integer :: Nb
contains

    subroutine x0_y0_z0()
        implicit none
        integer :: k
        real :: x1,x2,x3
        call random_seed()
        do k=1,N,1
            call random_number(x1)
            call random_number(x2)
            call random_number(x3)
            x(k)=0.001*(k-N/2.0)+x1
            y(k)=0.002*(k-N/2.0)+x2
            z(k)=0.003*(k-N/2.0)+x3
            write(101,*) x(k),y(k),z(k)
        end do
    end subroutine x0_y0_z0

    subroutine fnf(p)
        implicit none
        integer :: i,j,p,t
        real :: func_x(N),func_y(N),func_z(N),func_e(N),a,b,d,k1,k2,k3,k4,II,JJ(N),C(N),r,s,x_i0,E_ext,x_s,&
                &lambda,thet_s,C1(N),C2(N),Li(N),Em,f
        k1=0.7
        k2=0.001
        k3=0.0
        k4=9.0
        a=1.0
        b=3.0
        d=5.0
        r=0.01
        s=5.0
        x_s=2.0
        x_i0=-1.6
        lambda=10.0
        thet_s=-0.25
        E_ext=0.0
        E=0.0
        Em=1.5
        f=12.0
        II=3.5
        do Nb=1,N,1
            do t=1,MaxT,1
                JJ=0.0
                C1=0.0
                C2=0.0
                C=0.0
                do i=1,N,1
                    do j=1,N,1
                        JJ(i)=JJ(i)+neighbour_matrix_2(i,j)*(x(j)-x(i))
                        C1(i)=C1(i)+neighbour_matrix_2P(i,j)*(1.0/(1.0+exp(-lambda*(x(j)-thet_s))))
                        C2(i)=C2(i)+neighbour_matrix_2(i,j)*(1.0/(1.0+exp(-lambda*(x(j)-thet_s))))
                    end do
                    C(i)=C1(i)-C2(i)
                    func_x(i)=x(i)+h*(y(i)-a*x(i)**3+b*x(i)**2-z(i)+II+k3*JJ(i)+((-k4/(2.0*p-2.0))*(x_s-x(i))*C(i)))
                    func_y(i)=y(i)+h*(1.0-d*x(i)**2-y(i)+k1*E(i))
                    func_z(i)=z(i)+h*(r*(s*(x(i)-x_i0)-z(i)))
                    if(i>=Nb) then
                        func_e(i)=E(i)+h*(k2*y(i)+Em*sin(2.0*PI*f*t*h))
                    else
                        func_e(i)=0.0
                    end if
                    if(Nb==11) then
                        write(107,*) t,x(i)
                    end if
                    xx(i,t)=x(i)
                end do
                x=func_x
                y=func_y
                z=func_z
                E=func_e
            end do
            call SIDM(xx)
        end do
        return
    end subroutine fnf

    !SI不相干强度
    subroutine SIDM(xx)
        implicit none
        integer :: i,j,t,m
        real(kind=8) :: delta
        real(kind=8) :: xx(N,MaxT),sig(Nb),S(Nb),ZZm(N)
        real(kind=8) :: Sm,SI,DM,ZZ,nn
        delta=0.02*abs(maxval(xx(:,:))-minval(xx(:,:)))
        nn=real(N)/real(Nb)
        SI=0.0
        DM=0.0
        ZZm=0.0
        do i=1,N-1,1
            do t=MaxT/2,MaxT,1
                ZZm(i)=ZZm(i)+(xx(i,t)-xx(i+1,t))
            end do
        end do
        do t=MaxT/2,MaxT,1
            ZZm(N)=ZZm(N)+(xx(N,t)-xx(1,t))
        end do
        ZZm=ZZm/real(MaxT/2)
        sig=0.0
        do m=1,Nb,1
            do t=MaxT/2,MaxT,1
                ZZ=0.0
                do j=nn*(m-1)+1,m*nn,1
                    if(j/=N) then
                        ZZ=ZZ+((xx(j,t)-xx(j+1,t))-ZZm(m))**2
                    else
                        ZZ=ZZ+((xx(N,t)-xx(1,t))-ZZm(m))**2
                    end if
                end do
                sig(m)=sig(m)+sqrt((1.0/nn)*ZZ)
            end do
        end do
        sig=sig/real(MaxT/2)
        do m=1,Nb,1
            Sm=delta-sig(m)
            write(*,*) sig(m),delta
            if(Sm>0.0) then
                S(m)=1.0
            else
                S(m)=0.0
            end if
        end do
        SI=1.0-(1.0/Nb)*sum(S(:))
        do m=1,Nb-1,1
            DM=DM+abs(S(m)-S(m+1))
        end do
        if(Nb==N) then
            DM=DM+abs(S(N)-S(1))
        end if
        DM=DM/2.0
        write(105,*) Nb,SI
        write(106,*) Nb,DM
        write(*,*) Nb,SI,DM,sum(S(:))
        return
    end subroutine SIDM

    !建立网状结构

    !k=2
    subroutine neighbour_2(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix_2(N,N))
        !初始化网络矩阵
        neighbour_matrix_2=0
        do i=1,N-1,1
            do j=i+1,N,1
                if(abs(i-j)<=K/2.or.abs(i-j)>=N-K/2) then
                    neighbour_matrix_2(i,j)=1
                    neighbour_matrix_2(j,i)=1
                end if
            end do
        end do
        !输出矩阵
        do i=1,N,1
            do j=1,N,1
                write(103,"(I2)",advance='no') neighbour_matrix_2(i,j)
            end do
            write(103,*)
        end do
        return
    end subroutine neighbour_2

    !k=2*p+2
    subroutine neighbour_2P(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix_2P(N,N))
        !初始化网络矩阵
        neighbour_matrix_2P=0
        do i=1,N-1,1
            do j=i+1,N,1
                if((abs(i-j)<=K/2.or.abs(i-j)>=N-K/2).and.abs(i-j)/=1) then
                    neighbour_matrix_2P(i,j)=1
                    neighbour_matrix_2P(j,i)=1
                end if
            end do
        end do
        !输出矩阵
        do i=1,N,1
            do j=1,N,1
                write(104,"(I2)",advance='no') neighbour_matrix_2P(i,j)
            end do
            write(104,*)
        end do
        return
    end subroutine neighbour_2P

end module HR



program main
    use HR
    implicit none
    integer :: p
    p=8
    open(101,file="data_x0_y0_z0.txt")
    open(102,file="data_i_t_x.txt")
    open(103,file="neighbour_matrix_2.txt")
    open(104,file="neighbour_matrix_2P.txt")
    open(105,file="SI.txt")
    open(106,file="DM.txt")
    open(107,file="Nb_11.txt")
    call x0_y0_z0()
    call neighbour_2(2)
    call neighbour_2P(2*p+2)
    call fnf(p)
    deallocate(neighbour_matrix_2)
    deallocate(neighbour_matrix_2P)
    close(101)
    close(102)
    close(103)
    close(104)
    close(105)
    close(106)
end program main

figure 8_A

figure 8_B

figure 8_C

figure 8_D

figure 8_E

figure 8_F

figure 8_G

figure 8_H

figure 9

module HR
    implicit none
    real,parameter :: h=0.01,PI=3.1415926
    integer,parameter :: N=100,MaxT=660000
    integer,allocatable :: neighbour_matrix_2(:,:),neighbour_matrix_2P(:,:)
    real :: x(N),y(N),z(N),E(N)
contains

    subroutine x0_y0_z0()
        implicit none
        integer :: k
        real :: x1,x2,x3
        call random_seed()
        do k=1,N,1
            call random_number(x1)
            call random_number(x2)
            call random_number(x3)
            x(k)=0.001*(k-N/2.0)+x1
            y(k)=0.002*(k-N/2.0)+x2
            z(k)=0.003*(k-N/2.0)+x3
            read(101,*) x(k),y(k),z(k)
        end do
    end subroutine x0_y0_z0

    subroutine fnf(p)
        implicit none
        integer :: i,j,p,t
        real :: func_x(N),func_y(N),func_z(N),func_e(N),a,b,d,k1,k2,k3,k4,II,JJ(N),C(N),r,s,x_i0,E_ext,x_s,&
                &lambda,thet_s,C1(N),C2(N),Li(N),Em,f
        k1=0.7
        k2=0.001
        k3=0.0
        k4=9.0
        a=1.0
        b=3.0
        d=5.0
        r=0.01
        s=5.0
        x_s=2.0
        x_i0=-1.6
        lambda=10.0
        thet_s=-0.25
        E_ext=0.0
        E=0.0
        Em=1.5
        f=12.0
        II=35.0
        do t=1,MaxT,1
            if(mod(t,5000)==0) then
                write(*,*) t
            end if
            JJ=0.0
            C1=0.0
            C2=0.0
            C=0.0
            do i=1,N,1
                do j=1,N,1
                    JJ(i)=JJ(i)+neighbour_matrix_2(i,j)*(x(j)-x(i))
                    C1(i)=C1(i)+neighbour_matrix_2P(i,j)*(1.0/(1.0+exp(-lambda*(x(j)-thet_s))))
                    C2(i)=C2(i)+neighbour_matrix_2(i,j)*(1.0/(1.0+exp(-lambda*(x(j)-thet_s))))
                end do
                C(i)=C1(i)-C2(i)
                func_x(i)=x(i)+h*(y(i)-a*x(i)**3+b*x(i)**2-z(i)+II+k3*JJ(i)+((-k4/(2.0*p-2.0))*(x_s-x(i))*C(i)))
                func_y(i)=y(i)+h*(1.0-d*x(i)**2-y(i)+k1*E(i))
                func_z(i)=z(i)+h*(r*(s*(x(i)-x_i0)-z(i)))
                if(i>50) then
                    func_e(i)=E(i)+h*(k2*y(i)+Em*sin(2.0*PI*f*t*h))

                else
                    func_e(i)=0.0
                end if
                if(t>600000.and.mod(t,50)==0) then
                    write(102,*) i,t-600000,x(i)
                end if
                if(t>600000.and.i==5) then
                    write(105,*) t-600000,x(i)
                end if   
                if(t>600000.and.i==55) then
                    write(106,*) t-600000,x(i)
                end if
            end do
            x=func_x
            y=func_y
            z=func_z
            E=func_e
        end do
        return
    end subroutine fnf

    !建立网状结构

    !k=2
    subroutine neighbour_2(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix_2(N,N))
        !初始化网络矩阵
        neighbour_matrix_2=0
        do i=1,N-1,1
            do j=i+1,N,1
                if(abs(i-j)<=K/2.or.abs(i-j)>=N-K/2) then
                    neighbour_matrix_2(i,j)=1
                    neighbour_matrix_2(j,i)=1
                end if
            end do
        end do
        !输出矩阵
        do i=1,N,1
            do j=1,N,1
                write(103,"(I2)",advance='no') neighbour_matrix_2(i,j)
            end do
            write(103,*)
        end do
        return
    end subroutine neighbour_2

    !k=2*p+2
    subroutine neighbour_2P(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix_2P(N,N))
        !初始化网络矩阵
        neighbour_matrix_2P=0
        do i=1,N-1,1
            do j=i+1,N,1
                if((abs(i-j)<=K/2.or.abs(i-j)>=N-K/2).and.abs(i-j)/=1) then
                    neighbour_matrix_2P(i,j)=1
                    neighbour_matrix_2P(j,i)=1
                end if
            end do
        end do
        !输出矩阵
        do i=1,N,1
            do j=1,N,1
                write(104,"(I2)",advance='no') neighbour_matrix_2P(i,j)
            end do
            write(104,*)
        end do
        return
    end subroutine neighbour_2P

end module HR



program main
    use HR
    implicit none
    integer :: p
    p=40
    open(101,file="data_x0_y0_z0.txt")
    open(102,file="data_i_t_x.txt")
    open(103,file="neighbour_matrix_2.txt")
    open(104,file="neighbour_matrix_2P.txt")
    open(105,file="t_x_5.txt")
    open(106,file="t_x_55.txt")
    call x0_y0_z0()
    call neighbour_2(2)
    call neighbour_2P(2*p+2)
    call fnf(p)
    deallocate(neighbour_matrix_2)
    deallocate(neighbour_matrix_2P)
    close(101)
    close(102)
    close(103)
    close(104)
end program main

figure 9_A

figure 9_B

figure 10

module HR
    implicit none
    real,parameter :: h=0.01,PI=3.1415926
    integer,parameter :: N=100,MaxT=160000
    integer,allocatable :: neighbour_matrix_2(:,:),neighbour_matrix_2P(:,:)
    real :: x(N),y(N),z(N),E(N)
contains

    subroutine x0_y0_z0()
        implicit none
        integer :: k
        real :: x1,x2,x3
        call random_seed()
        do k=1,N,1
            call random_number(x1)
            call random_number(x2)
            call random_number(x3)
            x(k)=0.001*(k-N/2.0)+x1
            y(k)=0.002*(k-N/2.0)+x2
            z(k)=0.003*(k-N/2.0)+x3
            read(101,*) x(k),y(k),z(k)
        end do
    end subroutine x0_y0_z0

    subroutine fnf(p)
        implicit none
        integer :: i,j,p,t
        real :: func_x(N),func_y(N),func_z(N),func_e(N),a,b,d,k1,k2,k3,k4,II,JJ(N),C(N),r,s,x_i0,E_ext,x_s,&
                &lambda,thet_s,C1(N),C2(N),Li(N),Em,f
        k1=0.7
        k2=0.001
        k3=2.0
        k4=3.0
        a=1.0
        b=3.0
        d=5.0
        r=0.01
        s=5.0
        x_s=2.0
        x_i0=-1.6
        lambda=10.0
        thet_s=-0.25
        E_ext=0.0
        E=0.0
        Em=1.5
        f=12.0
        II=3.5
        do t=1,MaxT,1
            if(mod(t,5000)==0) then
                write(*,*) t
            end if
            JJ=0.0
            C1=0.0
            C2=0.0
            C=0.0
            do i=1,N,1
                do j=1,N,1
                    JJ(i)=JJ(i)+neighbour_matrix_2(i,j)*(x(j)-x(i))
                    C1(i)=C1(i)+neighbour_matrix_2P(i,j)*(1.0/(1.0+exp(-lambda*(x(j)-thet_s))))
                    C2(i)=C2(i)+neighbour_matrix_2(i,j)*(1.0/(1.0+exp(-lambda*(x(j)-thet_s))))
                end do
                C(i)=C1(i)-C2(i)
                func_x(i)=x(i)+h*(y(i)-a*x(i)**3+b*x(i)**2-z(i)+II+k3*JJ(i)+((-k4/(2.0*p-2.0))*(x_s-x(i))*C(i)))
                func_y(i)=y(i)+h*(1.0-d*x(i)**2-y(i)+k1*E(i))
                func_z(i)=z(i)+h*(r*(s*(x(i)-x_i0)-z(i)))
                if(i<25.or.(i>50.and.i<75)) then
                    func_e(i)=E(i)+h*(k2*y(i)+Em*sin(2.0*PI*f*t*h))
                else
                    func_e(i)=0.0
                end if
                if(t>100000.and.mod(t,200)==0) then
                    write(102,*) i,t-100000,x(i)
                end if
                if(t>100000.and.i==10) then
                    write(105,*) t-100000,x(i)
                end if   
                if(t>100000.and.i==85) then
                    write(106,*) t-100000,x(i)
                end if
            end do
            x=func_x
            y=func_y
            z=func_z
            E=func_e
        end do
        return
    end subroutine fnf

    !建立网状结构

    !k=2
    subroutine neighbour_2(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix_2(N,N))
        !初始化网络矩阵
        neighbour_matrix_2=0
        do i=1,N-1,1
            do j=i+1,N,1
                if(abs(i-j)<=K/2.or.abs(i-j)>=N-K/2) then
                    neighbour_matrix_2(i,j)=1
                    neighbour_matrix_2(j,i)=1
                end if
            end do
        end do
        !输出矩阵
        do i=1,N,1
            do j=1,N,1
                write(103,"(I2)",advance='no') neighbour_matrix_2(i,j)
            end do
            write(103,*)
        end do
        return
    end subroutine neighbour_2

    !k=2*p+2
    subroutine neighbour_2P(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix_2P(N,N))
        !初始化网络矩阵
        neighbour_matrix_2P=0
        do i=1,N-1,1
            do j=i+1,N,1
                if((abs(i-j)<=K/2.or.abs(i-j)>=N-K/2).and.abs(i-j)/=1) then
                    neighbour_matrix_2P(i,j)=1
                    neighbour_matrix_2P(j,i)=1
                end if
            end do
        end do
        !输出矩阵
        do i=1,N,1
            do j=1,N,1
                write(104,"(I2)",advance='no') neighbour_matrix_2P(i,j)
            end do
            write(104,*)
        end do
        return
    end subroutine neighbour_2P

end module HR



program main
    use HR
    implicit none
    integer :: p
    p=40
    open(101,file="data_x0_y0_z0.txt")
    open(102,file="data_i_t_x.txt")
    open(103,file="neighbour_matrix_2.txt")
    open(104,file="neighbour_matrix_2P.txt")
    open(105,file="t_x_5.txt")
    open(106,file="t_x_55.txt")
    call x0_y0_z0()
    call neighbour_2(2)
    call neighbour_2P(2*p+2)
    call fnf(p)
    deallocate(neighbour_matrix_2)
    deallocate(neighbour_matrix_2P)
    close(101)
    close(102)
    close(103)
    close(104)
end program main

figure 10_A

figure 10_B

figure 11

• k1=0.7,k2=0.001,f=12

module HR
    implicit none
    real,parameter :: h=0.01,PI=3.1415926
    integer,parameter :: N=100,MaxT=390000
    integer,allocatable :: neighbour_matrix_2(:,:),neighbour_matrix_2P(:,:)
    real :: x(N),y(N),z(N),E(N)
contains

    subroutine x0_y0_z0()
        implicit none
        integer :: k
        real :: x1,x2,x3
        call random_seed()
        do k=1,N,1
            call random_number(x1)
            call random_number(x2)
            call random_number(x3)
            x(k)=0.001*(k-N/2.0)+x1/100.0
            y(k)=0.002*(k-N/2.0)+x2/100.0
            z(k)=0.003*(k-N/2.0)+x3/100.0
            write(101,*) x(k),y(k),z(k)
        end do
    end subroutine x0_y0_z0

    subroutine fnf(p)
        implicit none
        integer :: i,j,p,t
        real :: func_x(N),func_y(N),func_z(N),func_e(N),a,b,d,k1,k2,k3,k4,II,JJ(N),C(N),r,s,x_i0,E_ext,x_s,&
                &lambda,thet_s,C1(N),C2(N),Li(N),Em,f
        k1=0.7
        k2=0.001
        k3=0.0
        k4=9.0
        a=1.0
        b=3.0
        d=5.0
        r=0.01
        s=5.0
        x_s=2.0
        x_i0=-1.6
        lambda=10.0
        thet_s=-0.25
        E_ext=0.0
        E=0.0
        Em=1.5
        f=12.0
        II=3.5
        do t=1,MaxT,1
            if(mod(t,5000)==0) then
                write(*,*) t
            end if
            JJ=0.0
            C1=0.0
            C2=0.0
            C=0.0
            do i=1,N,1
                do j=1,N,1
                    JJ(i)=JJ(i)+neighbour_matrix_2(i,j)*(x(j)-x(i))
                    C1(i)=C1(i)+neighbour_matrix_2P(i,j)*(1.0/(1.0+exp(-lambda*(x(j)-thet_s))))
                    C2(i)=C2(i)+neighbour_matrix_2(i,j)*(1.0/(1.0+exp(-lambda*(x(j)-thet_s))))
                end do
                C(i)=C1(i)-C2(i)
                func_x(i)=x(i)+h*(y(i)-a*x(i)**3+b*x(i)**2-z(i)+II+k3*JJ(i)+((-k4/(2.0*p-2.0))*(x_s-x(i))*C(i)))
                func_y(i)=y(i)+h*(1.0-d*x(i)**2-y(i)+k1*E(i))
                func_z(i)=z(i)+h*(r*(s*(x(i)-x_i0)-z(i)))
                E_ext=Em*sin(2.0*PI*f*(t))
                func_e(i)=k2*y(i)+E_ext
                if(t>300000.and.mod(t,1000)==0) then
                    write(102,*) i,t-300000,x(i)
                end if  
            end do
            x=func_x
            y=func_y
            z=func_z
            E=func_e
        end do
        return
    end subroutine fnf

    !建立网状结构

    !k=2
    subroutine neighbour_2(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix_2(N,N))
        !初始化网络矩阵
        neighbour_matrix_2=0
        do i=1,N-1,1
            do j=i+1,N,1
                if(abs(i-j)<=K/2.or.abs(i-j)>=N-K/2) then
                    neighbour_matrix_2(i,j)=1
                    neighbour_matrix_2(j,i)=1
                end if
            end do
        end do
        !输出矩阵
        do i=1,N,1
            do j=1,N,1
                write(103,"(I2)",advance='no') neighbour_matrix_2(i,j)
            end do
            write(103,*)
        end do
        return
    end subroutine neighbour_2

    !k=2*p+2
    subroutine neighbour_2P(K)
        implicit none
        integer :: K !单节点K个邻居
        integer :: i,j,number
        allocate(neighbour_matrix_2P(N,N))
        !初始化网络矩阵
        neighbour_matrix_2P=0
        do i=1,N-1,1
            do j=i+1,N,1
                if((abs(i-j)<=K/2.or.abs(i-j)>=N-K/2).and.abs(i-j)/=1) then
                    neighbour_matrix_2P(i,j)=1
                    neighbour_matrix_2P(j,i)=1
                end if
            end do
        end do
        !输出矩阵
        do i=1,N,1
            do j=1,N,1
                write(104,"(I2)",advance='no') neighbour_matrix_2P(i,j)
            end do
            write(104,*)
        end do
        return
    end subroutine neighbour_2P

end module HR



program main
    use HR
    implicit none
    integer :: p
    p=2
    open(101,file="data_x0_y0_z0.txt")
    open(102,file="data_i_t_x.txt")
    open(103,file="neighbour_matrix_2.txt")
    open(104,file="neighbour_matrix_2P.txt")
    call x0_y0_z0()
    call neighbour_2(2)
    call neighbour_2P(2*p+2)
    call fnf(p)
    deallocate(neighbour_matrix_2)
    deallocate(neighbour_matrix_2P)
    close(101)
    close(102)
    close(103)
    close(104)
end program main

figure 11_A

figure 11_B

figure 11_C

figure 11_D

figure 11_E

figure 11_F

figure 12

figure 13